Main Page   Class Hierarchy   Compound List   File List   Compound Members   File Members  

SparseMat.cc

Go to the documentation of this file.
00001 /*
00002     File:           SparseMat.cc
00003 
00004     Function:       Implements SparseMat.h
00005 
00006     Author(s):      Andrew Willmott
00007 
00008     Copyright:      (c) 1995-2000, Andrew Willmott
00009 
00010     Notes:          
00011 
00012 */
00013 
00014 
00015 #include "vl/SparseMat.h"
00016 #include <iomanip.h>
00017 #include <ctype.h>
00018 #include <string.h>
00019 #include <stdarg.h>
00020 #include "cl/Array.h"
00021 
00022 
00023 // --- SparseMat Constructors & Destructors -----------------------------------
00024 
00025 
00026 TSparseMat::TSparseMat() : rows(0), cols(0), row(0)
00027 {
00028 }
00029 
00030 TSparseMat::TSparseMat(Int nrows, Int ncols) : row(0)
00031 {
00032     SetSize(nrows, ncols);
00033     MakeZero(); // Adds in sentinels.
00034 }
00035 
00036 TSparseMat::TSparseMat(Int rows, Int cols, ZeroOrOne k) : row(0)
00037 {
00038     SetSize(rows, cols);
00039     MakeDiag(k);
00040 }
00041 
00042 TSparseMat::TSparseMat(Int rows, Int cols, Block k) : row(0)
00043 {
00044     SetSize(rows, cols);
00045     MakeBlock(k);
00046 }
00047 
00048 TSparseMat::TSparseMat(const TSparseMat &m) : row(0)
00049 {
00050     Assert(m.row != 0, "(SparseMat) Can't construct from null matrix");
00051     
00052     Int i;
00053     
00054     SetSize(m.rows, m.cols);
00055     
00056     for (i = 0; i < rows; i++)
00057         row[i] = m.row[i];
00058 }
00059 
00060 TSparseMat::TSparseMat(const TSubSMat &m) : row(0)
00061 {
00062     SELF = m;
00063 }
00064 
00065 TSparseMat::TSparseMat(const TMat &m) : row(0)
00066 {
00067     Int i;
00068     
00069     SetSize(m.Rows(), m.Cols());
00070 
00071     for (i = 0; i < rows; i++)
00072         row[i] = m[i];
00073 }
00074 
00075 TSparseMat::~TSparseMat()
00076 {
00077     delete[] row;
00078 }
00079 
00080 Void TSparseMat::SetSize(Int m, Int n)
00081 {
00082     Int i;
00083 
00084     Assert(m > 0 && n > 0, "(SparseMat::SetSize) illegal matrix size");
00085 
00086     rows = m;
00087     cols = n;
00088     
00089     delete[] row;   
00090     row = new TMSparseVec[rows];
00091     
00092     Assert(row != 0, "(SparseMat::SetSize) Out of memory");
00093     
00094     for (i=0; i < rows; i++)
00095         row[i].SetSize(cols);
00096 }
00097 
00098 
00099 // --- SparseMat Assignment Operators -----------------------------------------
00100 
00101 
00102 TSparseMat &TSparseMat::operator = (const TSparseMat &m)
00103 {   
00104     Int i;
00105     
00106     if (rows == 0)
00107         SetSize(m.Rows(), m.Cols());
00108     else
00109         Assert(rows == m.Rows() && cols == m.Cols(),
00110              "(SparseMat::=) matrices not the same size");
00111         
00112     for (i = 0; i < rows; i++)
00113         row[i] = m.row[i];
00114     
00115     return(SELF);
00116 }
00117       
00118 TSparseMat &TSparseMat::operator = (const TMat &m)
00119 {   
00120     Int i;
00121     
00122     SetSize(m.Rows(), m.Cols());
00123 
00124     for (i = 0; i < rows; i++)
00125         row[i] = m[i];
00126 
00127     return(SELF);
00128 }
00129 
00130 TSparseMat &TSparseMat::operator = (const TSubSMat &m)
00131 {   
00132     Int i;
00133     
00134     SetSize(m.Rows(), m.Cols());
00135 
00136     for (i = 0; i < rows; i++)
00137         row[i] = m[i];
00138 
00139     return(SELF);
00140 }
00141 
00142 Void TSparseMat::MakeZero()
00143 {
00144     Int     i;
00145     
00146     for (i = 0; i < rows; i++)
00147         row[i] = vl_zero;       
00148 }
00149 
00150 Void TSparseMat::MakeDiag(TMReal k)
00151 {
00152     Int     i;
00153     
00154     for (i = 0; i < rows; i++)
00155         row[i].MakeUnit(i, k);      
00156 }
00157 
00158 Void TSparseMat::MakeBlock(TMReal k)
00159 {
00160     Int     i;
00161     
00162     for (i = 0; i < rows; i++)
00163         row[i].MakeBlock(k);
00164 }
00165 
00166 
00167 // --- Mat Assignment Operators -----------------------------------------------
00168 
00169 
00170 TSparseMat &operator += (TSparseMat &m, const TSparseMat &n)
00171 {
00172     Assert(n.Rows() == m.Rows(), "(SparseMat::+=) matrix rows don't match");    
00173     
00174     Int     i;
00175     
00176     for (i = 0; i < m.Rows(); i++) 
00177         m[i] += n[i];
00178     
00179     return(m);
00180 }
00181 
00182 TSparseMat &operator -= (TSparseMat &m, const TSparseMat &n)
00183 {
00184     Assert(n.Rows() == m.Rows(), "(SparseMat::-=) matrix rows don't match");    
00185     
00186     Int     i;
00187     
00188     for (i = 0; i < m.Rows(); i++) 
00189         m[i] -= n[i];
00190     
00191     return(m);
00192 }
00193 
00194 TSparseMat &operator *= (TSparseMat &m, const TSparseMat &n)
00195 {
00196     Assert(m.Cols() == n.Cols(), "(SparseMat::*=) matrix columns don't match"); 
00197     
00198     Int     i;
00199     
00200     for (i = 0; i < m.Rows(); i++) 
00201         m[i] *= (TSparseMat &) n;
00202     
00203     return(m);
00204 }
00205 
00206 TSparseMat &operator *= (TSparseMat &m, TMReal s)
00207 {   
00208     Int     i;
00209     
00210     for (i = 0; i < m.Rows(); i++) 
00211         m[i] *= s;
00212     
00213     return(m);
00214 }
00215 
00216 TSparseMat &operator /= (TSparseMat &m, TMReal s)
00217 {   
00218     Int     i;
00219     
00220     for (i = 0; i < m.Rows(); i++) 
00221         m[i] /= s;
00222     
00223     return(m);
00224 }
00225 
00226 // --- Mat Comparison Operators -----------------------------------------------
00227 
00228 Bool operator == (const TSparseMat &m, const TSparseMat &n)
00229 {
00230     Assert(n.Rows() == m.Rows(), "(SparseMat::==) matrix rows don't match");    
00231     
00232     Int     i;
00233     
00234     for (i = 0; i < m.Rows(); i++) 
00235         if (m[i] != n[i])
00236             return(0);
00237 
00238     return(1);
00239 }
00240 
00241 Bool operator != (const TSparseMat &m, const TSparseMat &n)
00242 {
00243     Assert(n.Rows() == m.Rows(), "(SparseMat::!=) matrix rows don't match");    
00244     
00245     Int     i;
00246     
00247     for (i = 0; i < m.Rows(); i++) 
00248         if (m[i] != n[i])
00249             return(1);
00250 
00251     return(0);
00252 }
00253 
00254 // --- Mat Arithmetic Operators -----------------------------------------------
00255 
00256 TSparseMat operator + (const TSparseMat &m, const TSparseMat &n)
00257 {
00258     Assert(n.Rows() == m.Rows(), "(SparseMat::+) matrix rows don't match"); 
00259     
00260     TSparseMat  result(m.Rows(), m.Cols());
00261     Int     i;
00262     
00263     for (i = 0; i < m.Rows(); i++) 
00264         result[i] = m[i] + n[i];
00265     
00266     return(result);
00267 }
00268 
00269 TSparseMat operator - (const TSparseMat &m, const TSparseMat &n)
00270 {
00271     Assert(n.Rows() == m.Rows(), "(SparseMat::-) matrix rows don't match"); 
00272     
00273     TSparseMat  result(m.Rows(), m.Cols());
00274     Int     i;
00275     
00276     for (i = 0; i < m.Rows(); i++) 
00277         result[i] = m[i] - n[i];
00278     
00279     return(result);
00280 }
00281 
00282 TSparseMat operator - (const TSparseMat &m)
00283 {
00284     TSparseMat  result(m.Rows(), m.Cols());
00285     Int     i;
00286     
00287     for (i = 0; i < m.Rows(); i++) 
00288         result[i] = -m[i];
00289     
00290     return(result);
00291 }
00292 
00293 TSparseMat operator * (const TSparseMat &m, const TSparseMat &n)
00294 {
00295     Assert(m.Cols() == n.Rows(), "(SparseMat::*m) matrix cols don't match");    
00296     
00297     TSparseMat  result(m.Rows(), n.Cols());
00298     Int     i;
00299     
00300     for (i = 0; i < m.Rows(); i++) 
00301         result[i] = m[i] * n;
00302     
00303     return(result);
00304 }
00305 
00306 TSparseVec operator * (const TSparseMat &m, const TSparseVec &v)
00307 {
00308     Assert(m.Cols() == v.Elts(), "(SparseMat::m*v) matrix and vector sizes don't match");
00309     
00310     Int         i;
00311     TSparseVec  result(m.Rows());
00312     
00313     for (i = 0; i < m.Rows(); i++) 
00314         result.AddElt(i, dot(m[i], v));
00315     
00316     result.End();
00317         
00318     return(result);
00319 }
00320 
00321 TMVec operator * (const TSparseMat &m, const TMVec &v)
00322 {
00323     Assert(m.Cols() == v.Elts(), "(SparseMat::*v) matrix and vector sizes don't match");
00324     
00325     Int         i;
00326     TMVec       result(m.Rows());
00327     
00328     for (i = 0; i < m.Rows(); i++) 
00329         result[i] = dot(m[i], v);
00330     
00331     return(result);
00332 }
00333 
00334 TSparseMat operator * (const TSparseMat &m, TMReal s)
00335 {
00336     Int     i;
00337     TSparseMat  result(m.Rows(), m.Cols());
00338     
00339     for (i = 0; i < m.Rows(); i++) 
00340         result[i] = m[i] * s;
00341     
00342     return(result);
00343 }
00344 
00345 TSparseMat operator / (const TSparseMat &m, TMReal s)
00346 {
00347     Int     i;
00348     TSparseMat  result(m.Rows(), m.Cols());
00349     
00350     for (i = 0; i < m.Rows(); i++) 
00351         result[i] = m[i] / s;
00352     
00353     return(result);
00354 }
00355 
00356 
00357 // --- Mat-Vec Functions ------------------------------------------------------
00358 
00359 
00360 TMSparseVec operator * (const TSparseVec &v, const TSparseMat &m)           // v * m
00361 {
00362     Assert(v.Elts() == m.Rows(), "(Mat::v*m) vector/matrix sizes don't match");
00363     
00364     TMSparseVec result(m.Cols());
00365     Int         i;
00366     TSVIter     j(v);
00367 
00368     result = vl_zero;   
00369     
00370     // v * M = v[0] * m[0] + v[1] * m[1] + ...
00371     
00372     for (i = 0; i < m.Rows(); i++) 
00373     {       
00374         j.Inc(i);
00375         if (j.Exists(i))
00376             result += m[i] * j.Data();
00377     }
00378     
00379     return(result);
00380 }
00381 
00382 TMVec operator * (const TMVec &v, const TSparseMat &m)          // v * m
00383 {
00384     Assert(v.Elts() == m.Rows(), "(Mat::v*m) vector/matrix sizes don't match");
00385     
00386     TMVec       result(m.Cols());
00387     Int         i, j;
00388 
00389     result = vl_zero;   
00390 
00391     // v * M = v[0] * m[0] + v[1] * m[1] + ...
00392     
00393     for (i = 0; i < m.Rows(); i++) 
00394         if (len(v[i]) > 0)
00395             for (j = 0; j < m[i].NumItems() - 1; j++)
00396                 result[m[i][j].index] += v[i] * m[i][j].elt;
00397     
00398     return(result);
00399 }
00400 
00401 TSparseVec &operator *= (TSparseVec &v, const TSparseMat &m)        // v *= m
00402 {
00403     TSparseVec t = v * m;
00404     v.Replace(t);
00405     return(v);
00406 }
00407 
00408 TMVec &operator *= (TMVec &v, const TSparseMat &m)                  // v *= m
00409 {
00410     v = v * m;
00411     return(v);
00412 }
00413 
00414 
00415 // --- Mat Special Functions --------------------------------------------------
00416 
00417 
00418 TSparseMat trans(const TSparseMat &m)
00419 {
00420     Int         i, j;
00421     TSparseMat  result(m.Cols(), m.Rows());
00422         
00423     for (i = 0; i < result.Rows(); i++)
00424         result[i].Begin();
00425 
00426     for (i = 0; i < m.Rows(); i++)
00427     {
00428         // For row i of 'm', add its elements to the ends of the
00429         // appropriate row of 'result'.
00430     
00431         for (j = 0; j < m[i].NumItems() - 1; j++)
00432             result[m[i][j].index].AddNZElt(i, m[i][j].elt);
00433     }
00434     
00435     for (i = 0; i < result.Rows(); i++)
00436         result[i].End();
00437             
00438     return(result);
00439 }
00440 
00441 TMReal trace(const TSparseMat &m)
00442 {
00443     Int     i;
00444     TSVIter j;
00445     TMReal  result = vl_0;
00446     
00447     for (i = 0; i < m.Rows(); i++) 
00448     {
00449         j.Begin(m[i]);
00450         
00451         // Find element i of m[i], and if it exists,
00452         // add it to the result.
00453         
00454         if (j.IncTo(i))
00455             result += j.Data();
00456     }
00457     
00458     return(result);
00459 }
00460 
00461 TSparseMat oprod(const TSparseVec &a, const TSparseVec &b)
00462 // returns outerproduct of a and b:  a * trans(b)
00463 {
00464     TSparseMat  result;
00465     TSVIter     i;
00466     
00467     result.SetSize(a.Elts(), b.Elts());
00468     result = vl_0;
00469         
00470     for (i.Begin(a); !i.AtEnd(); i.Inc())
00471     {
00472         result[i.Index()] = b;
00473         result[i.Index()] *= i.Data();
00474     }
00475     
00476     return(result);
00477 }
00478 
00479 TSparseMat oprods(const TVec &a, const TVec &b)
00480 // returns outerproduct of a and b:  a * trans(b)
00481 {
00482     TSparseMat  result;
00483     Int         i;
00484 
00485     result.SetSize(a.Elts(), b.Elts());
00486     
00487     for (i = 0; i < a.Elts(); i++)
00488         if (TSparseVec::IsNonZero(a[i]))
00489         {
00490             result[i] = b;
00491             result[i] *= a[i];
00492         }
00493         else
00494             result[i] = vl_0;
00495 
00496     return(result);
00497 }
00498 
00499 
00500 // --- Mat Input & Output -----------------------------------------------------
00501 
00502 
00503 ostream &operator << (ostream &s, const TSparseMat &m)
00504 {
00505     Int i, w = s.width();
00506 
00507     s << '[';
00508     for (i = 0; i < m.Rows() - 1; i++)
00509         s << setw(w) << m[i] << endl;
00510     s << setw(w) << m[i] << ']' << endl;
00511     
00512     return(s);
00513 }
00514 
00515 istream &operator >> (istream &s, TSparseMat &m)
00516 {
00517     Array< TMSparseVec > array; 
00518     Int     i;
00519     
00520     s >> array;                     // Read input into array of SparseVecs
00521     
00522     m.SetSize(array.NumItems(), array[0].NumItems());
00523     
00524     for (i = 0; i < m.Rows(); i++)  // copy the result into m
00525     {
00526         Assert(m.Cols() == array[i].NumItems(), "(Mat/>>) different sized matrix rows");
00527         m[i] = array[i];
00528     }
00529     
00530     return(s);
00531 }
00532 
00533 TSparseMat inv(const TSparseMat &m, TMReal *determinant, TMReal pEps)
00534 {
00535     Assert(m.IsSquare(), "(inv) Matrix not square"); 
00536 
00537     //  Note that the sparse version of inv() is actually simpler than
00538     //  the normal version.
00539 
00540     Int             i, j, k;
00541     Int             n = m.Rows();
00542     TMReal          t, det, pivot;
00543     Real            max;
00544     TSparseMat      A(m);
00545     TSparseMat      B(n, n, vl_I);      
00546 
00547     // ---------- Forward elimination ---------- ------------------------------
00548     
00549     det = vl_1;
00550     
00551     for (i = 0; i < n; i++)     
00552     {           
00553         // Eliminate in column i, below diagonal
00554         
00555         max = -1.0;
00556         
00557         // Find a pivot for column i 
00558         // For the sparse rows, we take advantage of the fact that if A(i, k) exists,
00559         // it will be the first element in the sparse vector. (Previous elements will have
00560         // been zeroed by previous elimination steps.)
00561         
00562         for (k = i; k < n; k++) 
00563             if ((A[k][0].index == i) && (len(A[k][0].elt) > max))
00564             {
00565                 pivot = A[k][0].elt;
00566                 max = len(pivot);
00567                 j = k;
00568             }
00569             
00570         // If no nonzero pivot exists, A must be singular...
00571 
00572         Assert(max > pEps, "(inv) Matrix is singular");
00573             
00574         // Swap rows i and j
00575 
00576         if (j != i)
00577         {   
00578             //  Sparse matrix rows are just arrays: we can swap the contents
00579             //  of the two arrays efficiently by using the array Swap function.
00580 
00581             A[i].SwapWith(A[j]);
00582             B[i].SwapWith(B[j]);
00583                 
00584             det = -det;
00585         }
00586         
00587         det *= pivot;
00588         
00589         A[i] /= pivot;
00590         B[i] /= pivot;    
00591            
00592         for (j = i + 1; j < n; j++)
00593         {
00594             // Eliminate in rows below i 
00595             // Again, if A[j,i] exists, it will be the first non-zero element of the row.
00596             
00597             if (A[j][0].index == i)
00598             {
00599                 t = A[j][0].elt;
00600                 A[j] -= A[i] * t;
00601                 B[j] -= B[i] * t;
00602             }
00603         }
00604     }
00605 
00606     // ---------- Backward elimination ---------- -----------------------------
00607 
00608     for (i = 1; i < n; i++)         // Eliminate in column i, above diag 
00609     {       
00610         for (j = 0; j < i; j++)         // Eliminate in rows above i 
00611         {       
00612             if (A[j][1].index == i)
00613             {
00614                 t = A[j][1].elt;
00615                 A[j] -= A[i] * t;
00616                 B[j] -= B[i] * t;
00617             }
00618         }
00619     }
00620     if (determinant)
00621         *determinant = det;
00622     return(B);
00623 }
00624 

Generated at Sat Aug 5 00:16:47 2000 for Class Library by doxygen 1.1.0 written by Dimitri van Heesch, © 1997-2000