NEWMAT7.CPP Source File Main Page Namespace List Class Hierarchy Alphabetical List Compound List File List Namespace Members Compound Members File Members
NEWMAT7.CPP

Go to the documentation of this file.
00001 //Owner: Fred
00002 //$$ newmat7.cpp     Invert, solve, binary operations
00003 
00004 // Copyright (C) 1991,2,3,4: R B Davies
00005 
00006 #include "include.h"
00007 
00008 #include "newmat.h"
00009 #include "newmatrc.h"
00010 
00011 #ifdef use_namespace
00012 namespace NEWMAT {
00013 #endif
00014 
00015 #ifdef DO_REPORT
00016 #define REPORT { static ExeCounter ExeCount(__LINE__,7); ++ExeCount; }
00017 #else
00018 #define REPORT {}
00019 #endif
00020 
00021     /***************************** solve routines ******************************/
00022 
00023     GeneralMatrix* GeneralMatrix::MakeSolver() {
00024         REPORT
00025             GeneralMatrix* gm = new CroutMatrix(*this);
00026         MatrixErrorNoSpace(gm); gm->ReleaseAndDelete(); return gm;
00027     }
00028 
00029     GeneralMatrix* Matrix::MakeSolver() {
00030         REPORT
00031             GeneralMatrix* gm = new CroutMatrix(*this);
00032         MatrixErrorNoSpace(gm); gm->ReleaseAndDelete(); return gm;
00033     }
00034 
00035     void CroutMatrix::Solver(MatrixColX& mcout, const MatrixColX& mcin) {
00036         REPORT
00037             int i = mcin.skip; Real* el = mcin.data-i; Real* el1 = el;
00038             while (i--) *el++ = 0.0;
00039             el += mcin.storage; i = nrows - mcin.skip - mcin.storage;
00040             while (i--) *el++ = 0.0;
00041             lubksb(el1, mcout.skip);
00042     }
00043 
00044     // Do we need check for entirely zero output?
00045 
00046     void UpperTriangularMatrix::Solver(MatrixColX& mcout,
00047                                        const MatrixColX& mcin) {
00048         REPORT
00049             int i = mcin.skip-mcout.skip; Real* elx = mcin.data-i;
00050         while (i-- > 0) *elx++ = 0.0;
00051         int nr = mcin.skip+mcin.storage;
00052         elx = mcin.data+mcin.storage; Real* el = elx;
00053         int j = mcout.skip+mcout.storage-nr; int nc = ncols-nr; i = nr-mcout.skip;
00054         while (j-- > 0) *elx++ = 0.0;
00055         Real* Ael = store + (nr*(2*ncols-nr+1))/2; j = 0;
00056         while (i-- > 0) {
00057             elx = el; Real sum = 0.0; int jx = j++; Ael -= nc;
00058             while (jx--) sum += *(--Ael) * *(--elx);
00059             elx--; *elx = (*elx - sum) / *(--Ael);
00060         }
00061     }
00062 
00063     void LowerTriangularMatrix::Solver(MatrixColX& mcout,
00064                                        const MatrixColX& mcin) {
00065         REPORT
00066             int i = mcin.skip-mcout.skip; Real* elx = mcin.data-i;
00067         while (i-- > 0) *elx++ = 0.0;
00068         int nc = mcin.skip; i = nc+mcin.storage; elx = mcin.data+mcin.storage;
00069         int nr = mcout.skip+mcout.storage; int j = nr-i; i = nr-nc;
00070         while (j-- > 0) *elx++ = 0.0;
00071         Real* el = mcin.data; Real* Ael = store + (nc*(nc+1))/2; j = 0;
00072         while (i-- > 0) {
00073             elx = el; Real sum = 0.0; int jx = j++; Ael += nc;
00074             while (jx--) sum += *Ael++ * *elx++;
00075             *elx = (*elx - sum) / *Ael++;
00076         }
00077     }
00078 
00079     /******************* carry out binary operations *************************/
00080 
00081     static GeneralMatrix*
00082     GeneralAdd(GeneralMatrix*,GeneralMatrix*,AddedMatrix*,MatrixType);
00083     static GeneralMatrix*
00084     GeneralSub(GeneralMatrix*,GeneralMatrix*,SubtractedMatrix*,MatrixType);
00085     static GeneralMatrix*
00086     GeneralMult(GeneralMatrix*,GeneralMatrix*,MultipliedMatrix*,MatrixType);
00087     static GeneralMatrix*
00088     GeneralSolv(GeneralMatrix*,GeneralMatrix*,BaseMatrix*,MatrixType);
00089     static GeneralMatrix*
00090     GeneralSP(GeneralMatrix*,GeneralMatrix*,SPMatrix*,MatrixType);
00091     static GeneralMatrix*
00092     GeneralElmDivide(GeneralMatrix*,GeneralMatrix*,ElmDivideMatrix*,MatrixType);
00093 
00094     GeneralMatrix* AddedMatrix::Evaluate(MatrixType mt) {
00095         REPORT
00096             gm1=((BaseMatrix*&)bm1)->Evaluate();
00097         gm2=((BaseMatrix*&)bm2)->Evaluate();
00098 #ifdef TEMPS_DESTROYED_QUICKLY
00099         GeneralMatrix* gmx;
00100         Try { gmx = GeneralAdd(gm1,gm2,this,mt); }
00101         CatchAll { delete this; ReThrow; }
00102         delete this; return gmx;
00103 #else
00104         return GeneralAdd(gm1,gm2,this,mt);
00105 #endif
00106     }
00107 
00108     GeneralMatrix* SubtractedMatrix::Evaluate(MatrixType mt) {
00109         REPORT
00110             gm1=((BaseMatrix*&)bm1)->Evaluate();
00111         gm2=((BaseMatrix*&)bm2)->Evaluate();
00112 #ifdef TEMPS_DESTROYED_QUICKLY
00113         GeneralMatrix* gmx;
00114         Try { gmx = GeneralSub(gm1,gm2,this,mt); }
00115         CatchAll { delete this; ReThrow; }
00116         delete this; return gmx;
00117 #else
00118         return GeneralSub(gm1,gm2,this,mt);
00119 #endif
00120     }
00121 
00122     GeneralMatrix* MultipliedMatrix::Evaluate(MatrixType mt) {
00123         REPORT
00124             gm2 = ((BaseMatrix*&)bm2)->Evaluate();
00125         gm2 = gm2->Evaluate(gm2->Type().MultRHS());   // no symmetric on RHS
00126         gm1=((BaseMatrix*&)bm1)->Evaluate();
00127 #ifdef TEMPS_DESTROYED_QUICKLY
00128         GeneralMatrix* gmx;
00129         Try { gmx = GeneralMult(gm1, gm2, this, mt); }
00130         CatchAll { delete this; ReThrow; }
00131         delete this; return gmx;
00132 #else
00133         return GeneralMult(gm1, gm2, this, mt);
00134 #endif
00135     }
00136 
00137     GeneralMatrix* SolvedMatrix::Evaluate(MatrixType mt) {
00138         REPORT
00139             gm1=((BaseMatrix*&)bm1)->Evaluate();
00140         gm2=((BaseMatrix*&)bm2)->Evaluate();
00141 #ifdef TEMPS_DESTROYED_QUICKLY
00142         GeneralMatrix* gmx;
00143         Try { gmx = GeneralSolv(gm1,gm2,this,mt); }
00144         CatchAll { delete this; ReThrow; }
00145         delete this; return gmx;
00146 #else
00147         return GeneralSolv(gm1,gm2,this,mt);
00148 #endif
00149     }
00150 
00151     GeneralMatrix* SPMatrix::Evaluate(MatrixType mt) {
00152         REPORT
00153             gm1=((BaseMatrix*&)bm1)->Evaluate();
00154         gm2=((BaseMatrix*&)bm2)->Evaluate();
00155 #ifdef TEMPS_DESTROYED_QUICKLY
00156         GeneralMatrix* gmx;
00157         Try { gmx = GeneralSP(gm1,gm2,this,mt); }
00158         CatchAll { delete this; ReThrow; }
00159         delete this; return gmx;
00160 #else
00161         return GeneralSP(gm1,gm2,this,mt);
00162 #endif
00163     }
00164 
00165     GeneralMatrix* ElmDivideMatrix::Evaluate(MatrixType mt) {
00166         REPORT
00167             gm1=((BaseMatrix*&)bm1)->Evaluate();
00168         gm2=((BaseMatrix*&)bm2)->Evaluate();
00169 #ifdef TEMPS_DESTROYED_QUICKLY
00170         GeneralMatrix* gmx;
00171         Try { gmx = GeneralElmDivide(gm1,gm2,this,mt); }
00172         CatchAll { delete this; ReThrow; }
00173         delete this; return gmx;
00174 #else
00175         return GeneralElmDivide(gm1,gm2,this,mt);
00176 #endif
00177     }
00178 
00179     // routines for adding or subtracting matrices of identical storage structure
00180 
00181     static void Add(GeneralMatrix* gm, GeneralMatrix* gm1, GeneralMatrix* gm2) {
00182         REPORT
00183             Real* s1=gm1->Store(); Real* s2=gm2->Store();
00184         Real* s=gm->Store(); int i=gm->Storage() >> 2;
00185         while (i--) {
00186             *s++ = *s1++ + *s2++; *s++ = *s1++ + *s2++;
00187             *s++ = *s1++ + *s2++; *s++ = *s1++ + *s2++;
00188         }
00189         i=gm->Storage() & 3; while (i--) *s++ = *s1++ + *s2++;
00190     }
00191 
00192     static void Add(GeneralMatrix* gm, GeneralMatrix* gm2) {
00193         REPORT
00194             Real* s2=gm2->Store(); Real* s=gm->Store(); int i=gm->Storage() >> 2;
00195             while (i--)
00196             { *s++ += *s2++; *s++ += *s2++; *s++ += *s2++; *s++ += *s2++; }
00197             i=gm->Storage() & 3; while (i--) *s++ += *s2++;
00198     }
00199 
00200     static void Subtract(GeneralMatrix* gm, GeneralMatrix* gm1, GeneralMatrix* gm2) {
00201         REPORT
00202             Real* s1=gm1->Store(); Real* s2=gm2->Store();
00203         Real* s=gm->Store(); int i=gm->Storage() >> 2;
00204         while (i--) {
00205             *s++ = *s1++ - *s2++; *s++ = *s1++ - *s2++;
00206             *s++ = *s1++ - *s2++; *s++ = *s1++ - *s2++;
00207         }
00208         i=gm->Storage() & 3; while (i--) *s++ = *s1++ - *s2++;
00209     }
00210 
00211     static void Subtract(GeneralMatrix* gm, GeneralMatrix* gm2) {
00212         REPORT
00213             Real* s2=gm2->Store(); Real* s=gm->Store(); int i=gm->Storage() >> 2;
00214             while (i--)
00215             { *s++ -= *s2++; *s++ -= *s2++; *s++ -= *s2++; *s++ -= *s2++; }
00216             i=gm->Storage() & 3; while (i--) *s++ -= *s2++;
00217     }
00218 
00219     static void ReverseSubtract(GeneralMatrix* gm, GeneralMatrix* gm2) {
00220         REPORT
00221             Real* s2=gm2->Store(); Real* s=gm->Store(); int i=gm->Storage() >> 2;
00222             while (i--) {
00223                 *s = *s2++ - *s; s++; *s = *s2++ - *s; s++;
00224                 *s = *s2++ - *s; s++; *s = *s2++ - *s; s++;
00225             }
00226             i=gm->Storage() & 3; while (i--) { *s = *s2++ - *s; s++; }
00227     }
00228 
00229     static void SP(GeneralMatrix* gm, GeneralMatrix* gm1, GeneralMatrix* gm2) {
00230         REPORT
00231             Real* s1=gm1->Store(); Real* s2=gm2->Store();
00232         Real* s=gm->Store(); int i=gm->Storage() >> 2;
00233         while (i--) {
00234             *s++ = *s1++ * *s2++; *s++ = *s1++ * *s2++;
00235             *s++ = *s1++ * *s2++; *s++ = *s1++ * *s2++;
00236         }
00237         i=gm->Storage() & 3; while (i--) *s++ = *s1++ * *s2++;
00238     }
00239 
00240     static void SP(GeneralMatrix* gm, GeneralMatrix* gm2) {
00241         REPORT
00242             Real* s2=gm2->Store(); Real* s=gm->Store(); int i=gm->Storage() >> 2;
00243             while (i--)
00244             { *s++ *= *s2++; *s++ *= *s2++; *s++ *= *s2++; *s++ *= *s2++; }
00245             i=gm->Storage() & 3; while (i--) *s++ *= *s2++;
00246     }
00247 
00248     static void ElmDivide(GeneralMatrix* gm, GeneralMatrix* gm1, GeneralMatrix* gm2) {
00249         REPORT
00250             Real* s1=gm1->Store(); Real* s2=gm2->Store();
00251         Real* s=gm->Store(); int i=gm->Storage() >> 2;
00252         while (i--) {
00253             *s++ = *s1++ / *s2++; *s++ = *s1++ / *s2++;
00254             *s++ = *s1++ / *s2++; *s++ = *s1++ / *s2++;
00255         }
00256         i=gm->Storage() & 3; while (i--) *s++ = *s1++ / *s2++;
00257     }
00258 
00259     static void ElmDivide(GeneralMatrix* gm, GeneralMatrix* gm2) {
00260         REPORT
00261             Real* s2=gm2->Store(); Real* s=gm->Store(); int i=gm->Storage() >> 2;
00262             while (i--)
00263             { *s++ /= *s2++; *s++ /= *s2++; *s++ /= *s2++; *s++ /= *s2++; }
00264             i=gm->Storage() & 3; while (i--) *s++ /= *s2++;
00265     }
00266 
00267     // routines for adding or subtracting matrices of different storage structure
00268 
00269     static void AddDS(GeneralMatrix* gm, GeneralMatrix* gm1, GeneralMatrix* gm2) {
00270         int nr = gm->Nrows();
00271         MatrixRow mr1(gm1, LoadOnEntry); MatrixRow mr2(gm2, LoadOnEntry);
00272         MatrixRow mr(gm, StoreOnExit+DirectPart);
00273         while (nr--) { mr.Add(mr1,mr2); mr1.Next(); mr2.Next(); mr.Next(); }
00274     }
00275 
00276     static void AddDS(GeneralMatrix* gm, GeneralMatrix* gm2) {
00277         // Add into first argument
00278         int nr = gm->Nrows();
00279         MatrixRow mr(gm, StoreOnExit+LoadOnEntry+DirectPart);
00280         MatrixRow mr2(gm2, LoadOnEntry);
00281         while (nr--) { mr.Add(mr2); mr.Next(); mr2.Next(); }
00282     }
00283 
00284     static void SubtractDS
00285     (GeneralMatrix* gm, GeneralMatrix* gm1, GeneralMatrix* gm2) {
00286         int nr = gm->Nrows();
00287         MatrixRow mr1(gm1, LoadOnEntry); MatrixRow mr2(gm2, LoadOnEntry);
00288         MatrixRow mr(gm, StoreOnExit+DirectPart);
00289         while (nr--) { mr.Sub(mr1,mr2); mr1.Next(); mr2.Next(); mr.Next(); }
00290     }
00291 
00292     static void SubtractDS(GeneralMatrix* gm, GeneralMatrix* gm2) {
00293         int nr = gm->Nrows();
00294         MatrixRow mr(gm, LoadOnEntry+StoreOnExit+DirectPart);
00295         MatrixRow mr2(gm2, LoadOnEntry);
00296         while (nr--) { mr.Sub(mr2); mr.Next(); mr2.Next(); }
00297     }
00298 
00299     static void ReverseSubtractDS(GeneralMatrix* gm, GeneralMatrix* gm2) {
00300         int nr = gm->Nrows();
00301         MatrixRow mr(gm, LoadOnEntry+StoreOnExit+DirectPart);
00302         MatrixRow mr2(gm2, LoadOnEntry);
00303         while (nr--) { mr.RevSub(mr2); mr2.Next(); mr.Next(); }
00304     }
00305 
00306     static void SPDS(GeneralMatrix* gm, GeneralMatrix* gm1, GeneralMatrix* gm2) {
00307         int nr = gm->Nrows();
00308         MatrixRow mr1(gm1, LoadOnEntry); MatrixRow mr2(gm2, LoadOnEntry);
00309         MatrixRow mr(gm, StoreOnExit+DirectPart);
00310         while (nr--) { mr.Multiply(mr1,mr2); mr1.Next(); mr2.Next(); mr.Next(); }
00311     }
00312 
00313     static void SPDS(GeneralMatrix* gm, GeneralMatrix* gm2) {
00314         // SP into first argument
00315         int nr = gm->Nrows();
00316         MatrixRow mr(gm, StoreOnExit+LoadOnEntry+DirectPart);
00317         MatrixRow mr2(gm2, LoadOnEntry);
00318         while (nr--) { mr.Multiply(mr2); mr.Next(); mr2.Next(); }
00319     }
00320 
00321 #ifdef __GNUG__
00322     void AddedMatrix::SelectVersion
00323     (MatrixType mtx, int& c1, int& c2) const
00324 #else
00325         void AddedMatrix::SelectVersion
00326     (MatrixType mtx, bool& c1, bool& c2) const
00327 #endif
00328     {
00329         // for determining version of add and subtract routines
00330         // will need to modify if further matrix structures are introduced
00331         MatrixBandWidth bw1 = gm1->BandWidth();
00332         MatrixBandWidth bw2 = gm2->BandWidth();
00333         MatrixBandWidth bwx(-1); if (mtx.IsBand()) bwx = bw1 + bw2;
00334         c1 = (mtx == gm1->Type()) && (bwx == bw1);
00335         c2 = (mtx == gm2->Type()) && (bwx == bw2);
00336     }
00337 
00338     static GeneralMatrix* GeneralAdd(GeneralMatrix* gm1, GeneralMatrix* gm2,
00339                                      AddedMatrix* am, MatrixType mtx) {
00340         Tracer tr("GeneralAdd");
00341         int nr=gm1->Nrows(); int nc=gm1->Ncols();
00342         if (nr!=gm2->Nrows() || nc!=gm2->Ncols())
00343             Throw(IncompatibleDimensionsException(*gm1, *gm2));
00344         Compare(gm1->Type() + gm2->Type(),mtx);
00345 #ifdef __GNUG__
00346         int c1,c2; am->SelectVersion(mtx,c1,c2);
00347 #else
00348         bool c1,c2; am->SelectVersion(mtx,c1,c2);     // causes problems for g++
00349 #endif
00350         if (c1 && c2) {
00351             if (gm1->reuse()) { REPORT Add(gm1,gm2); gm2->tDelete(); return gm1; }
00352             else if (gm2->reuse()) { REPORT Add(gm2,gm1); return gm2; }
00353             else {
00354                 REPORT GeneralMatrix* gmx = gm1->Type().New(nr,nc,am);
00355                 gmx->ReleaseAndDelete(); Add(gmx,gm1,gm2); return gmx;
00356             }
00357         }
00358         else {
00359             if (c1 && gm1->reuse() )                    // must have type test first
00360             { REPORT AddDS(gm1,gm2); gm2->tDelete(); return gm1; }
00361             else if (c2 && gm2->reuse() )
00362             { REPORT AddDS(gm2,gm1); if (!c1) gm1->tDelete(); return gm2; }
00363             else {
00364                 REPORT
00365                     GeneralMatrix* gmx = mtx.New(nr,nc,am); AddDS(gmx,gm1,gm2);
00366                 if (!c1) gm1->tDelete(); if (!c2) gm2->tDelete();
00367                 gmx->ReleaseAndDelete(); return gmx;
00368             }
00369         }
00370     }
00371 
00372     static GeneralMatrix* GeneralSub(GeneralMatrix* gm1, GeneralMatrix* gm2,
00373                                      SubtractedMatrix* sm, MatrixType mtx) {
00374         Tracer tr("GeneralSub");
00375         Compare(gm1->Type() + gm2->Type(),mtx);
00376         int nr=gm1->Nrows(); int nc=gm1->Ncols();
00377         if (nr!=gm2->Nrows() || nc!=gm2->Ncols())
00378             Throw(IncompatibleDimensionsException(*gm1, *gm2));
00379 #ifdef __GNUG__
00380         int c1,c2; sm->SelectVersion(mtx,c1,c2);
00381 #else
00382         bool c1,c2; sm->SelectVersion(mtx,c1,c2);     // causes problems for g++
00383 #endif
00384         if (c1 && c2) {
00385             if (gm1->reuse())
00386             { REPORT Subtract(gm1,gm2); gm2->tDelete(); return gm1; }
00387             else if (gm2->reuse()) { REPORT ReverseSubtract(gm2,gm1); return gm2; }
00388             else {
00389                 REPORT
00390                     GeneralMatrix* gmx = gm1->Type().New(nr,nc,sm);
00391                 gmx->ReleaseAndDelete(); Subtract(gmx,gm1,gm2); return gmx;
00392             }
00393         }
00394         else {
00395             if ( c1 && gm1->reuse() )
00396             { REPORT  SubtractDS(gm1,gm2); gm2->tDelete(); return gm1; }
00397             else if ( c2 && gm2->reuse() ) {
00398                 REPORT
00399                     ReverseSubtractDS(gm2,gm1); if (!c1) gm1->tDelete(); return gm2;
00400             }
00401             else {
00402                 REPORT
00403                     GeneralMatrix* gmx = mtx.New(nr,nc,sm); SubtractDS(gmx,gm1,gm2);
00404                 if (!c1) gm1->tDelete(); if (!c2) gm2->tDelete();
00405                 gmx->ReleaseAndDelete(); return gmx;
00406             }
00407         }
00408     }
00409 
00410     static GeneralMatrix* GeneralMult1(GeneralMatrix* gm1, GeneralMatrix* gm2,
00411                                        MultipliedMatrix* mm, MatrixType mtx) {
00412         REPORT
00413             Tracer tr("GeneralMult1");
00414         int nr=gm1->Nrows(); int nc=gm2->Ncols();
00415         if (gm1->Ncols() !=gm2->Nrows())
00416             Throw(IncompatibleDimensionsException(*gm1, *gm2));
00417         GeneralMatrix* gmx = mtx.New(nr,nc,mm);
00418 
00419         MatrixCol mcx(gmx, StoreOnExit+DirectPart);
00420         MatrixCol mc2(gm2, LoadOnEntry);
00421         while (nc--) {
00422             MatrixRow mr1(gm1, LoadOnEntry, mcx.Skip());
00423             Real* el = mcx.Data();                      // pointer to an element
00424             int n = mcx.Storage();
00425             while (n--) { *(el++) = DotProd(mr1,mc2); mr1.Next(); }
00426             mc2.Next(); mcx.Next();
00427         }
00428         gmx->ReleaseAndDelete(); gm1->tDelete(); gm2->tDelete(); return gmx;
00429     }
00430 
00431     static GeneralMatrix* GeneralMult2(GeneralMatrix* gm1, GeneralMatrix* gm2,
00432                                        MultipliedMatrix* mm, MatrixType mtx) {
00433         // version that accesses by row only - not good for thin matrices
00434         // or column vectors in right hand term. Needs fixing
00435         REPORT
00436             Tracer tr("GeneralMult2");
00437         int nr=gm1->Nrows(); int nc=gm2->Ncols();
00438         if (gm1->Ncols() !=gm2->Nrows())
00439             Throw(IncompatibleDimensionsException(*gm1, *gm2));
00440         GeneralMatrix* gmx = mtx.New(nr,nc,mm);
00441 
00442         Real* el = gmx->Store(); int n = gmx->Storage();
00443         while (n--) *el++ = 0.0;
00444         MatrixRow mrx(gmx, LoadOnEntry+StoreOnExit+DirectPart);
00445         MatrixRow mr1(gm1, LoadOnEntry);
00446         while (nr--) {
00447             MatrixRow mr2(gm2, LoadOnEntry, mr1.Skip());
00448             el = mr1.Data();                            // pointer to an element
00449             n = mr1.Storage();
00450             while (n--) { mrx.AddScaled(mr2, *el++); mr2.Next(); }
00451             mr1.Next(); mrx.Next();
00452         }
00453         gmx->ReleaseAndDelete(); gm1->tDelete(); gm2->tDelete(); return gmx;
00454     }
00455 
00456     static GeneralMatrix* mmMult(GeneralMatrix* gm1, GeneralMatrix* gm2) {
00457         // matrix multiplication for type Matrix only
00458         REPORT
00459             Tracer tr("MatrixMult");
00460 
00461         int nr=gm1->Nrows(); int ncr=gm1->Ncols(); int nc=gm2->Ncols();
00462         if (ncr != gm2->Nrows()) Throw(IncompatibleDimensionsException(*gm1,*gm2));
00463 
00464         Matrix* gm = new Matrix(nr,nc); MatrixErrorNoSpace(gm);
00465 
00466         Real* s1=gm1->Store(); Real* s2=gm2->Store(); Real* s=gm->Store();
00467 
00468         if (ncr) {
00469             while (nr--) {
00470                 Real* s2x = s2; int j = ncr;
00471                 Real* sx = s; Real f = *s1++; int k = nc;
00472                 while (k--) *sx++ = f * *s2x++;
00473                 while (--j)
00474                 { sx = s; f = *s1++; k = nc; while (k--) *sx++ += f * *s2x++; }
00475                 s = sx;
00476             }
00477         }
00478         else *gm = 0.0;
00479 
00480         gm->ReleaseAndDelete(); gm1->tDelete(); gm2->tDelete(); return gm;
00481     }
00482 
00483     static GeneralMatrix* GeneralMult(GeneralMatrix* gm1, GeneralMatrix* gm2,
00484                                       MultipliedMatrix* mm, MatrixType mtx) {
00485         if ( Rectangular(gm1->Type(), gm2->Type(), mtx)) return mmMult(gm1, gm2);
00486         else {
00487             Compare(gm1->Type() * gm2->Type(),mtx);
00488             int nr = gm2->Nrows(); int nc = gm2->Ncols();
00489             if (nc <= 5 && nr > nc) return GeneralMult1(gm1, gm2, mm, mtx);
00490             else return GeneralMult2(gm1, gm2, mm, mtx);
00491         }
00492     }
00493 
00494 #ifdef __GNUG__
00495     void SPMatrix::SelectVersion
00496     (MatrixType mtx, int& c1, int& c2) const
00497 #else
00498         void SPMatrix::SelectVersion
00499     (MatrixType mtx, bool& c1, bool& c2) const
00500 #endif
00501     {
00502         // for determining version of SP routines
00503         // will need to modify if further matrix structures are introduced
00504         MatrixBandWidth bw1 = gm1->BandWidth();
00505         MatrixBandWidth bw2 = gm2->BandWidth();
00506         MatrixBandWidth bwx(-1);  if (mtx.IsBand()) bwx = bw1.minimum(bw2);
00507         c1 = (mtx == gm1->Type()) && (bwx == bw1);
00508         c2 = (mtx == gm2->Type()) && (bwx == bw2);
00509     }
00510 
00511     static GeneralMatrix* GeneralSP(GeneralMatrix* gm1, GeneralMatrix* gm2,
00512                                     SPMatrix* am, MatrixType mtx) {
00513         Tracer tr("GeneralSP");
00514         int nr=gm1->Nrows(); int nc=gm1->Ncols();
00515         if (nr!=gm2->Nrows() || nc!=gm2->Ncols())
00516             Throw(IncompatibleDimensionsException(*gm1, *gm2));
00517         Compare(gm1->Type().SP(gm2->Type()),mtx);
00518 #ifdef __GNUG__
00519         int c1,c2; am->SelectVersion(mtx,c1,c2);
00520 #else
00521         bool c1,c2; am->SelectVersion(mtx,c1,c2);     // causes problems for g++
00522 #endif
00523         if (c1 && c2) {
00524             if (gm1->reuse()) { REPORT SP(gm1,gm2); gm2->tDelete(); return gm1; }
00525             else if (gm2->reuse()) { REPORT SP(gm2,gm1); return gm2; }
00526             else {
00527                 REPORT GeneralMatrix* gmx = gm1->Type().New(nr,nc,am);
00528                 gmx->ReleaseAndDelete(); SP(gmx,gm1,gm2); return gmx;
00529             }
00530         }
00531         else {
00532             if (c1 && gm1->reuse() )                    // must have type test first
00533             { REPORT SPDS(gm1,gm2); gm2->tDelete(); return gm1; }
00534             else if (c2 && gm2->reuse() )
00535             { REPORT SPDS(gm2,gm1); if (!c1) gm1->tDelete(); return gm2; }
00536             else {
00537                 REPORT
00538                     GeneralMatrix* gmx = mtx.New(nr,nc,am); SPDS(gmx,gm1,gm2);
00539                 if (!c1) gm1->tDelete(); if (!c2) gm2->tDelete();
00540                 gmx->ReleaseAndDelete(); return gmx;
00541             }
00542         }
00543     }
00544 
00545     // fred 0512
00546     static GeneralMatrix* GeneralElmDivide(GeneralMatrix* gm1, GeneralMatrix* gm2,
00547                                            ElmDivideMatrix* am, MatrixType mtx) {
00548         Tracer tr("GeneralElmDivide");
00549         int nr=gm1->Nrows(); int nc=gm1->Ncols();
00550         if (nr!=gm2->Nrows() || nc!=gm2->Ncols())
00551             Throw(IncompatibleDimensionsException(*gm1, *gm2));
00552         Compare(gm1->Type().SP(gm2->Type()),mtx);     // BUG keep using SP's?
00553 #ifdef __GNUG__
00554         int c1,c2; am->SelectVersion(mtx,c1,c2);
00555 #else
00556         bool c1,c2; am->SelectVersion(mtx,c1,c2);     // causes problems for g++
00557 #endif
00558         if (c1 && c2) {
00559             if (gm1->reuse()) { REPORT ElmDivide(gm1,gm2); gm2->tDelete(); return gm1; }
00560             else if (gm2->reuse()) { REPORT ElmDivide(gm2,gm1); return gm2; }
00561             else {
00562                 REPORT GeneralMatrix* gmx = gm1->Type().New(nr,nc,am);
00563                 gmx->ReleaseAndDelete(); ElmDivide(gmx,gm1,gm2); return gmx;
00564             }
00565         }
00566         else {
00567             return gm1;                                 // BUG here
00568             /*
00569               if (c1 && gm1->reuse() )               // must have type test first
00570               {
00571               //REPORT SPDS(gm1,gm2); gm2->tDelete(); return gm1;
00572               }
00573               else if (c2 && gm2->reuse() )
00574               {
00575               //REPORT SPDS(gm2,gm1); if (!c1) gm1->tDelete(); return gm2;
00576               }
00577               else
00578               {
00579               REPORT
00580               GeneralMatrix* gmx = mtx.New(nr,nc,am); SPDS(gmx,gm1,gm2);
00581               if (!c1) gm1->tDelete(); if (!c2) gm2->tDelete();
00582               gmx->ReleaseAndDelete(); return gmx;
00583               }*/
00584         }
00585     }
00586 
00587     static GeneralMatrix* GeneralSolv(GeneralMatrix* gm1, GeneralMatrix* gm2,
00588                                       BaseMatrix* sm, MatrixType mtx) {
00589         REPORT
00590             Tracer tr("GeneralSolv");
00591         Compare(gm1->Type().i() * gm2->Type(),mtx);
00592         int nr=gm1->Nrows(); int nc=gm2->Ncols();
00593         if (gm1->Ncols() != gm2->Nrows())
00594             Throw(IncompatibleDimensionsException(*gm1, *gm2));
00595         GeneralMatrix* gmx = mtx.New(nr,nc,sm); MatrixErrorNoSpace(gmx);
00596         Real* r = new Real [nr]; MatrixErrorNoSpace(r);
00597 #ifndef ATandT
00598         MONITOR_REAL_NEW("Make   (GenSolv)",nr,r)
00599             // deleted for ATandT, to balance deletion below
00600 #endif
00601             GeneralMatrix* gms = gm1->MakeSolver();
00602         Try {
00603             // copy to and from r
00604             MatrixColX mcx(gmx, r, StoreOnExit+DirectPart);
00605             // this must be inside Try so mcx is destroyed before gmx
00606             MatrixColX mc2(gm2, r, LoadOnEntry);
00607             while (nc--) { gms->Solver(mcx, mc2); mcx.Next(); mc2.Next(); }
00608         }
00609         CatchAll {
00610             gms->tDelete();
00611             delete gmx;                                 // <--------------------
00612             gm2->tDelete();
00613 #ifndef ATandT
00614             MONITOR_REAL_DELETE("Delete (GenSolv)",nr,r)
00615                 // ATandT version 2.1 gives an internal error
00616 #endif
00617                 delete [] r;
00618             ReThrow;
00619         }
00620         gms->tDelete(); gmx->ReleaseAndDelete(); gm2->tDelete();
00621 #ifndef ATandT
00622         MONITOR_REAL_DELETE("Delete (GenSolv)",nr,r)
00623             // ATandT version 2.1 gives an internal error
00624 #endif
00625             delete [] r;
00626         return gmx;
00627     }
00628 
00629     GeneralMatrix* InvertedMatrix::Evaluate(MatrixType mtx) {
00630         // Matrix Inversion - use solve routines
00631         Tracer tr("InvertedMatrix::Evaluate");
00632         REPORT
00633             gm=((BaseMatrix*&)bm)->Evaluate();
00634         int n = gm->Nrows(); DiagonalMatrix I(n); I=1.0;
00635 #ifdef TEMPS_DESTROYED_QUICKLY
00636         GeneralMatrix* gmx;
00637         Try {
00638             Compare(gm->Type().i(),mtx);
00639             SkipConversionCheck SCC;                    // otherwise inverting
00640             // symmetric causes a problem
00641             gmx = GeneralSolv(gm,&I,this,mtx);
00642         }
00643         CatchAll { delete this; ReThrow; }
00644         delete this; return gmx;
00645 #else
00646         Compare(gm->Type().i(),mtx);
00647         SkipConversionCheck SCC;                      // otherwise inverting
00648         // symmetric causes a problem
00649         return GeneralSolv(gm,&I,this,mtx);
00650 #endif
00651     }
00652 
00653     /*************************** norm functions ****************************/
00654 
00655     Real BaseMatrix::Norm1() const {
00656         // maximum of sum of absolute values of a column
00657         REPORT
00658             GeneralMatrix* gm = ((BaseMatrix&)*this).Evaluate();
00659         int nc = gm->Ncols(); Real value = 0.0;
00660         MatrixCol mc(gm, LoadOnEntry);
00661         while (nc--)
00662         { Real v = mc.SumAbsoluteValue(); if (value < v) value = v; mc.Next(); }
00663         gm->tDelete(); return value;
00664     }
00665 
00666     Real BaseMatrix::NormInfinity() const {
00667         // maximum of sum of absolute values of a row
00668         REPORT
00669             GeneralMatrix* gm = ((BaseMatrix&)*this).Evaluate();
00670         int nr = gm->Nrows(); Real value = 0.0;
00671         MatrixRow mr(gm, LoadOnEntry);
00672         while (nr--)
00673         { Real v = mr.SumAbsoluteValue(); if (value < v) value = v; mr.Next(); }
00674         gm->tDelete(); return value;
00675     }
00676 
00677     /********************** Concatenation and stacking *************************/
00678 
00679     GeneralMatrix* ConcatenatedMatrix::Evaluate(MatrixType mtx) {
00680         REPORT
00681             Tracer tr("Concatenate");
00682 #ifdef TEMPS_DESTROYED_QUICKLY
00683         Try {
00684             gm2 = ((BaseMatrix*&)bm2)->Evaluate();
00685             gm1 = ((BaseMatrix*&)bm1)->Evaluate();
00686             Compare(gm1->Type() | gm2->Type(),mtx);
00687             int nr=gm1->Nrows(); int nc = gm1->Ncols() + gm2->Ncols();
00688             if (nr != gm2->Nrows())
00689                 Throw(IncompatibleDimensionsException(*gm1, *gm2));
00690             GeneralMatrix* gmx = mtx.New(nr,nc,this);
00691             MatrixRow mr1(gm1, LoadOnEntry); MatrixRow mr2(gm2, LoadOnEntry);
00692             MatrixRow mr(gmx, StoreOnExit+DirectPart);
00693             while (nr--) { mr.ConCat(mr1,mr2); mr1.Next(); mr2.Next(); mr.Next(); }
00694             gmx->ReleaseAndDelete(); gm1->tDelete(); gm2->tDelete(); delete this;
00695             return gmx;
00696         }
00697         CatchAll  { delete this; ReThrow; }
00698 #ifndef UseExceptions
00699         return 0;
00700 #endif
00701 #else
00702         gm2 = ((BaseMatrix*&)bm2)->Evaluate();
00703         gm1 = ((BaseMatrix*&)bm1)->Evaluate();
00704         Compare(gm1->Type() | gm2->Type(),mtx);
00705         int nr=gm1->Nrows(); int nc = gm1->Ncols() + gm2->Ncols();
00706         if (nr != gm2->Nrows())
00707             Throw(IncompatibleDimensionsException(*gm1, *gm2));
00708         GeneralMatrix* gmx = mtx.New(nr,nc,this);
00709         MatrixRow mr1(gm1, LoadOnEntry); MatrixRow mr2(gm2, LoadOnEntry);
00710         MatrixRow mr(gmx, StoreOnExit+DirectPart);
00711         while (nr--) { mr.ConCat(mr1,mr2); mr1.Next(); mr2.Next(); mr.Next(); }
00712         gmx->ReleaseAndDelete(); gm1->tDelete(); gm2->tDelete(); return gmx;
00713 #endif
00714     }
00715 
00716     GeneralMatrix* StackedMatrix::Evaluate(MatrixType mtx) {
00717         REPORT
00718             Tracer tr("Stack");
00719 #ifdef TEMPS_DESTROYED_QUICKLY
00720         Try {
00721             gm2 = ((BaseMatrix*&)bm2)->Evaluate();
00722             gm1 = ((BaseMatrix*&)bm1)->Evaluate();
00723             Compare(gm1->Type() & gm2->Type(),mtx);
00724             int nc=gm1->Ncols();
00725             int nr1 = gm1->Nrows(); int nr2 = gm2->Nrows();
00726             if (nc != gm2->Ncols())
00727                 Throw(IncompatibleDimensionsException(*gm1, *gm2));
00728             GeneralMatrix* gmx = mtx.New(nr1+nr2,nc,this);
00729             MatrixRow mr1(gm1, LoadOnEntry); MatrixRow mr2(gm2, LoadOnEntry);
00730             MatrixRow mr(gmx, StoreOnExit+DirectPart);
00731             while (nr1--) { mr.Copy(mr1); mr1.Next(); mr.Next(); }
00732             while (nr2--) { mr.Copy(mr2); mr2.Next(); mr.Next(); }
00733             gmx->ReleaseAndDelete(); gm1->tDelete(); gm2->tDelete(); delete this;
00734             return gmx;
00735         }
00736         CatchAll  { delete this; ReThrow; }
00737 #ifndef UseExceptions
00738         return 0;
00739 #endif
00740 #else
00741         gm2 = ((BaseMatrix*&)bm2)->Evaluate();
00742         gm1 = ((BaseMatrix*&)bm1)->Evaluate();
00743         Compare(gm1->Type() & gm2->Type(),mtx);
00744         int nc=gm1->Ncols();
00745         int nr1 = gm1->Nrows(); int nr2 = gm2->Nrows();
00746         if (nc != gm2->Ncols())
00747             Throw(IncompatibleDimensionsException(*gm1, *gm2));
00748         GeneralMatrix* gmx = mtx.New(nr1+nr2,nc,this);
00749         MatrixRow mr1(gm1, LoadOnEntry); MatrixRow mr2(gm2, LoadOnEntry);
00750         MatrixRow mr(gmx, StoreOnExit+DirectPart);
00751         while (nr1--) { mr.Copy(mr1); mr1.Next(); mr.Next(); }
00752         while (nr2--) { mr.Copy(mr2); mr2.Next(); mr.Next(); }
00753         gmx->ReleaseAndDelete(); gm1->tDelete(); gm2->tDelete(); return gmx;
00754 #endif
00755     }
00756 
00757     // ************************* equality of matrices ******************** //
00758 
00759     static bool RealEqual(Real* s1, Real* s2, int n) {
00760         int i = n >> 2;
00761         while (i--) {
00762             if (*s1++ != *s2++) return false; if (*s1++ != *s2++) return false;
00763             if (*s1++ != *s2++) return false; if (*s1++ != *s2++) return false;
00764         }
00765         i = n & 3; while (i--) if (*s1++ != *s2++) return false;
00766         return true;
00767     }
00768 
00769     static bool intEqual(int* s1, int* s2, int n) {
00770         int i = n >> 2;
00771         while (i--) {
00772             if (*s1++ != *s2++) return false; if (*s1++ != *s2++) return false;
00773             if (*s1++ != *s2++) return false; if (*s1++ != *s2++) return false;
00774         }
00775         i = n & 3; while (i--) if (*s1++ != *s2++) return false;
00776         return true;
00777     }
00778 
00779     bool operator==(const BaseMatrix& A, const BaseMatrix& B) {
00780         Tracer tr("BaseMatrix ==");
00781         REPORT
00782             GeneralMatrix* gmA = ((BaseMatrix&)A).Evaluate();
00783         GeneralMatrix* gmB = ((BaseMatrix&)B).Evaluate();
00784 
00785         if (gmA == gmB)                               // same matrix
00786         { REPORT gmA->tDelete(); return true; }
00787 
00788         if ( gmA->Nrows() != gmB->Nrows() || gmA->Ncols() != gmB->Ncols() )
00789             // different dimensions
00790         { REPORT gmA->tDelete(); gmB->tDelete(); return false; }
00791 
00792         // check for CroutMatrix or BandLUMatrix
00793         MatrixType AType = gmA->Type(); MatrixType BType = gmB->Type();
00794         if (AType.CannotConvert() || BType.CannotConvert() ) {
00795             REPORT
00796                 bool bx = gmA->IsEqual(*gmB);
00797             gmA->tDelete(); gmB->tDelete();
00798             return bx;
00799         }
00800 
00801         // is matrix storage the same
00802         // will need to modify if further matrix structures are introduced
00803         if (AType == BType && gmA->BandWidth() == gmB->BandWidth()) {
00804             // compare store
00805             REPORT
00806                 bool bx = RealEqual(gmA->Store(),gmB->Store(),gmA->Storage());
00807             gmA->tDelete(); gmB->tDelete();
00808             return bx;
00809         }
00810 
00811         // matrix storage different - just subtract
00812         REPORT  return IsZero(*gmA-*gmB);
00813     }
00814 
00815     bool operator==(const GeneralMatrix& A, const GeneralMatrix& B) {
00816         Tracer tr("GeneralMatrix ==");
00817         // May or may not call tDeletes
00818         REPORT
00819 
00820             if (&A == &B)                               // same matrix
00821             { REPORT return true; }
00822 
00823         if ( A.Nrows() != B.Nrows() || A.Ncols() != B.Ncols() ) {
00824             // different dimensions
00825             REPORT return false;
00826         }
00827 
00828         // check for CroutMatrix or BandLUMatrix
00829         MatrixType AType = A.Type(); MatrixType BType = B.Type();
00830         if (AType.CannotConvert() || BType.CannotConvert() )
00831         { REPORT  return A.IsEqual(B); }
00832 
00833         // is matrix storage the same
00834         // will need to modify if further matrix structures are introduced
00835         if (AType == BType && A.BandWidth() == B.BandWidth())
00836         { REPORT return RealEqual(A.Store(),B.Store(),A.Storage()); }
00837 
00838         // matrix storage different - just subtract
00839         REPORT  return IsZero(A-B);
00840     }
00841 
00842     bool GeneralMatrix::IsZero() const {
00843         REPORT
00844             Real* s=store; int i = storage >> 2;
00845         while (i--) {
00846             if (*s++) return false; if (*s++) return false;
00847             if (*s++) return false; if (*s++) return false;
00848         }
00849         i = storage & 3; while (i--) if (*s++) return false;
00850         return true;
00851     }
00852 
00853     bool IsZero(const BaseMatrix& A) {
00854         Tracer tr("BaseMatrix::IsZero");
00855         REPORT
00856             GeneralMatrix* gm1 = 0; bool bx;
00857         Try { gm1=((BaseMatrix&)A).Evaluate(); bx = gm1->IsZero(); }
00858         CatchAll { if (gm1) gm1->tDelete(); ReThrow; }
00859         gm1->tDelete();
00860         return bx;
00861     }
00862 
00863     // IsEqual functions - insist matrices are of same type
00864     // as well as equal values to be equal
00865 
00866     bool GeneralMatrix::IsEqual(const GeneralMatrix& A) const {
00867         Tracer tr("GeneralMatrix IsEqual");
00868         if (A.Type() != Type())                       // not same types
00869         { REPORT return false; }
00870         if (&A == this)                             // same matrix
00871         { REPORT  return true; }
00872         if (A.nrows != nrows || A.ncols != ncols)
00873             // different dimensions
00874         { REPORT return false; }
00875         // is matrix storage the same - compare store
00876         REPORT
00877             return RealEqual(A.store,store,storage);
00878     }
00879 
00880     bool CroutMatrix::IsEqual(const GeneralMatrix& A) const {
00881         Tracer tr("CroutMatrix IsEqual");
00882         if (A.Type() != Type())                       // not same types
00883         { REPORT return false; }
00884         if (&A == this)                             // same matrix
00885         { REPORT  return true; }
00886         if (A.nrows != nrows || A.ncols != ncols)
00887             // different dimensions
00888         { REPORT return false; }
00889         // is matrix storage the same - compare store
00890         REPORT
00891             return RealEqual(A.store,store,storage)
00892             && intEqual(((CroutMatrix&)A).indx, indx, nrows);
00893     }
00894 
00895     bool BandLUMatrix::IsEqual(const GeneralMatrix& A) const {
00896         Tracer tr("BandLUMatrix IsEqual");
00897         if (A.Type() != Type())                       // not same types
00898         { REPORT  return false; }
00899         if (&A == this)                             // same matrix
00900         { REPORT  return true; }
00901         if ( A.Nrows() != nrows || A.Ncols() != ncols
00902              || ((BandLUMatrix&)A).m1 != m1 || ((BandLUMatrix&)A).m2 != m2 )
00903             // different dimensions
00904         { REPORT  return false; }
00905 
00906         // matrix storage the same - compare store
00907         REPORT
00908             return RealEqual(A.Store(),store,storage)
00909             && RealEqual(((BandLUMatrix&)A).store2,store2,storage2)
00910             && intEqual(((BandLUMatrix&)A).indx, indx, nrows);
00911     }
00912 
00913 #ifdef use_namespace
00914 }
00915 #endif
Generated on Fri Aug 23 01:37:10 2002 for VirtualU by
1.2.17