newmat: newmat7.cpp Source File

00001 //$$ newmat7.cpp     Invert, solve, binary operations
00002 
00003 // Copyright (C) 1991,2,3,4: R B Davies
00004 
00005 #include "include.h"
00006 
00007 #include "newmat.h"
00008 #include "newmatrc.h"
00009 
00010 #ifdef use_namespace
00011 namespace NEWMAT {
00012 #endif
00013 
00014 
00015 #ifdef DO_REPORT
00016 #define REPORT { static ExeCounter ExeCount(__LINE__,7); ++ExeCount; }
00017 #else
00018 #define REPORT {}
00019 #endif
00020 
00021 
00022 //***************************** solve routines ******************************/
00023 
00024 GeneralMatrix* GeneralMatrix::MakeSolver()
00025 {
00026    REPORT
00027    GeneralMatrix* gm = new CroutMatrix(*this);
00028    MatrixErrorNoSpace(gm); gm->ReleaseAndDelete(); return gm;
00029 }
00030 
00031 GeneralMatrix* Matrix::MakeSolver()
00032 {
00033    REPORT
00034    GeneralMatrix* gm = new CroutMatrix(*this);
00035    MatrixErrorNoSpace(gm); gm->ReleaseAndDelete(); return gm;
00036 }
00037 
00038 void CroutMatrix::Solver(MatrixColX& mcout, const MatrixColX& mcin)
00039 {
00040    REPORT
00041    int i = mcin.skip; Real* el = mcin.data-i; Real* el1 = el;
00042    while (i--) *el++ = 0.0;
00043    el += mcin.storage; i = nrows_value - mcin.skip - mcin.storage;
00044    while (i--) *el++ = 0.0;
00045    lubksb(el1, mcout.skip);
00046 }
00047 
00048 
00049 // Do we need check for entirely zero output?
00050 
00051 void UpperTriangularMatrix::Solver(MatrixColX& mcout,
00052    const MatrixColX& mcin)
00053 {
00054    REPORT
00055    int i = mcin.skip-mcout.skip; Real* elx = mcin.data-i;
00056    while (i-- > 0) *elx++ = 0.0;
00057    int nr = mcin.skip+mcin.storage;
00058    elx = mcin.data+mcin.storage; Real* el = elx;
00059    int j = mcout.skip+mcout.storage-nr;
00060    int nc = ncols_value-nr; i = nr-mcout.skip;
00061    while (j-- > 0) *elx++ = 0.0;
00062    Real* Ael = store + (nr*(2*ncols_value-nr+1))/2; j = 0;
00063    while (i-- > 0)
00064    {
00065       elx = el; Real sum = 0.0; int jx = j++; Ael -= nc;
00066       while (jx--) sum += *(--Ael) * *(--elx);
00067       elx--; *elx = (*elx - sum) / *(--Ael);
00068    }
00069 }
00070 
00071 void LowerTriangularMatrix::Solver(MatrixColX& mcout,
00072    const MatrixColX& mcin)
00073 {
00074    REPORT
00075    int i = mcin.skip-mcout.skip; Real* elx = mcin.data-i;
00076    while (i-- > 0) *elx++ = 0.0;
00077    int nc = mcin.skip; i = nc+mcin.storage; elx = mcin.data+mcin.storage;
00078    int nr = mcout.skip+mcout.storage; int j = nr-i; i = nr-nc;
00079    while (j-- > 0) *elx++ = 0.0;
00080    Real* el = mcin.data; Real* Ael = store + (nc*(nc+1))/2; j = 0;
00081    while (i-- > 0)
00082    {
00083       elx = el; Real sum = 0.0; int jx = j++; Ael += nc;
00084       while (jx--) sum += *Ael++ * *elx++;
00085       *elx = (*elx - sum) / *Ael++;
00086    }
00087 }
00088 
00089 //******************* carry out binary operations *************************/
00090 
00091 static GeneralMatrix*
00092    GeneralMult(GeneralMatrix*,GeneralMatrix*,MultipliedMatrix*,MatrixType);
00093 static GeneralMatrix*
00094    GeneralSolv(GeneralMatrix*,GeneralMatrix*,BaseMatrix*,MatrixType);
00095 static GeneralMatrix*
00096    GeneralSolvI(GeneralMatrix*,BaseMatrix*,MatrixType);
00097 static GeneralMatrix*
00098    GeneralKP(GeneralMatrix*,GeneralMatrix*,KPMatrix*,MatrixType);
00099 
00100 GeneralMatrix* MultipliedMatrix::Evaluate(MatrixType mt)
00101 {
00102    REPORT
00103    gm2 = ((BaseMatrix*&)bm2)->Evaluate();
00104    gm2 = gm2->Evaluate(gm2->Type().MultRHS());     // no symmetric on RHS
00105    gm1=((BaseMatrix*&)bm1)->Evaluate();
00106    return GeneralMult(gm1, gm2, this, mt);
00107 }
00108 
00109 GeneralMatrix* SolvedMatrix::Evaluate(MatrixType mt)
00110 {
00111    REPORT
00112    gm1=((BaseMatrix*&)bm1)->Evaluate();
00113    gm2=((BaseMatrix*&)bm2)->Evaluate();
00114    return GeneralSolv(gm1,gm2,this,mt);
00115 }
00116 
00117 GeneralMatrix* KPMatrix::Evaluate(MatrixType mt)
00118 {
00119    REPORT
00120    gm1=((BaseMatrix*&)bm1)->Evaluate();
00121    gm2=((BaseMatrix*&)bm2)->Evaluate();
00122    return GeneralKP(gm1,gm2,this,mt);
00123 }
00124 
00125 // routines for adding or subtracting matrices of identical storage structure
00126 
00127 static void Add(GeneralMatrix* gm, GeneralMatrix* gm1, GeneralMatrix* gm2)
00128 {
00129    REPORT
00130    Real* s1=gm1->Store(); Real* s2=gm2->Store();
00131    Real* s=gm->Store(); int i=gm->Storage() >> 2;
00132    while (i--)
00133    {
00134        *s++ = *s1++ + *s2++; *s++ = *s1++ + *s2++;
00135        *s++ = *s1++ + *s2++; *s++ = *s1++ + *s2++;
00136    }
00137    i=gm->Storage() & 3; while (i--) *s++ = *s1++ + *s2++;
00138 }
00139 
00140 static void AddTo(GeneralMatrix* gm, const GeneralMatrix* gm2)
00141 {
00142    REPORT
00143    const Real* s2=gm2->Store(); Real* s=gm->Store(); int i=gm->Storage() >> 2;
00144    while (i--)
00145    { *s++ += *s2++; *s++ += *s2++; *s++ += *s2++; *s++ += *s2++; }
00146    i=gm->Storage() & 3; while (i--) *s++ += *s2++;
00147 }
00148 
00149 void GeneralMatrix::PlusEqual(const GeneralMatrix& gm)
00150 {
00151    REPORT
00152    if (nrows_value != gm.nrows_value || ncols_value != gm.ncols_value)
00153       Throw(IncompatibleDimensionsException(*this, gm));
00154    AddTo(this, &gm);
00155 }
00156 
00157 static void Subtract(GeneralMatrix* gm, GeneralMatrix* gm1, GeneralMatrix* gm2)
00158 {
00159    REPORT
00160    Real* s1=gm1->Store(); Real* s2=gm2->Store();
00161    Real* s=gm->Store(); int i=gm->Storage() >> 2;
00162    while (i--)
00163    {
00164        *s++ = *s1++ - *s2++; *s++ = *s1++ - *s2++;
00165        *s++ = *s1++ - *s2++; *s++ = *s1++ - *s2++;
00166    }
00167    i=gm->Storage() & 3; while (i--) *s++ = *s1++ - *s2++;
00168 }
00169 
00170 static void SubtractFrom(GeneralMatrix* gm, const GeneralMatrix* gm2)
00171 {
00172    REPORT
00173    const Real* s2=gm2->Store(); Real* s=gm->Store(); int i=gm->Storage() >> 2;
00174    while (i--)
00175    { *s++ -= *s2++; *s++ -= *s2++; *s++ -= *s2++; *s++ -= *s2++; }
00176    i=gm->Storage() & 3; while (i--) *s++ -= *s2++;
00177 }
00178 
00179 void GeneralMatrix::MinusEqual(const GeneralMatrix& gm)
00180 {
00181    REPORT
00182    if (nrows_value != gm.nrows_value || ncols_value != gm.ncols_value)
00183       Throw(IncompatibleDimensionsException(*this, gm));
00184    SubtractFrom(this, &gm);
00185 }
00186 
00187 static void ReverseSubtract(GeneralMatrix* gm, const GeneralMatrix* gm2)
00188 {
00189    REPORT
00190    const Real* s2=gm2->Store(); Real* s=gm->Store(); int i=gm->Storage() >> 2;
00191    while (i--)
00192    {
00193       *s = *s2++ - *s; s++; *s = *s2++ - *s; s++;
00194       *s = *s2++ - *s; s++; *s = *s2++ - *s; s++;
00195    }
00196    i=gm->Storage() & 3; while (i--) { *s = *s2++ - *s; s++; }
00197 }
00198 
00199 static void SP(GeneralMatrix* gm, GeneralMatrix* gm1, GeneralMatrix* gm2)
00200 {
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    {
00206        *s++ = *s1++ * *s2++; *s++ = *s1++ * *s2++;
00207        *s++ = *s1++ * *s2++; *s++ = *s1++ * *s2++;
00208    }
00209    i=gm->Storage() & 3; while (i--) *s++ = *s1++ * *s2++;
00210 }
00211 
00212 static void SP(GeneralMatrix* gm, GeneralMatrix* gm2)
00213 {
00214    REPORT
00215    Real* s2=gm2->Store(); Real* s=gm->Store(); int i=gm->Storage() >> 2;
00216    while (i--)
00217    { *s++ *= *s2++; *s++ *= *s2++; *s++ *= *s2++; *s++ *= *s2++; }
00218    i=gm->Storage() & 3; while (i--) *s++ *= *s2++;
00219 }
00220 
00221 // routines for adding or subtracting matrices of different storage structure
00222 
00223 static void AddDS(GeneralMatrix* gm, GeneralMatrix* gm1, GeneralMatrix* gm2)
00224 {
00225    REPORT
00226    int nr = gm->Nrows();
00227    MatrixRow mr1(gm1, LoadOnEntry); MatrixRow mr2(gm2, LoadOnEntry);
00228    MatrixRow mr(gm, StoreOnExit+DirectPart);
00229    while (nr--) { mr.Add(mr1,mr2); mr1.Next(); mr2.Next(); mr.Next(); }
00230 }
00231 
00232 static void AddDS(GeneralMatrix* gm, GeneralMatrix* gm2)
00233 // Add into first argument
00234 {
00235    REPORT
00236    int nr = gm->Nrows();
00237    MatrixRow mr(gm, StoreOnExit+LoadOnEntry+DirectPart);
00238    MatrixRow mr2(gm2, LoadOnEntry);
00239    while (nr--) { mr.Add(mr2); mr.Next(); mr2.Next(); }
00240 }
00241 
00242 static void SubtractDS
00243    (GeneralMatrix* gm, GeneralMatrix* gm1, GeneralMatrix* gm2)
00244 {
00245    REPORT
00246    int nr = gm->Nrows();
00247    MatrixRow mr1(gm1, LoadOnEntry); MatrixRow mr2(gm2, LoadOnEntry);
00248    MatrixRow mr(gm, StoreOnExit+DirectPart);
00249    while (nr--) { mr.Sub(mr1,mr2); mr1.Next(); mr2.Next(); mr.Next(); }
00250 }
00251 
00252 static void SubtractDS(GeneralMatrix* gm, GeneralMatrix* gm2)
00253 {
00254    REPORT
00255    int nr = gm->Nrows();
00256    MatrixRow mr(gm, LoadOnEntry+StoreOnExit+DirectPart);
00257    MatrixRow mr2(gm2, LoadOnEntry);
00258    while (nr--) { mr.Sub(mr2); mr.Next(); mr2.Next(); }
00259 }
00260 
00261 static void ReverseSubtractDS(GeneralMatrix* gm, GeneralMatrix* gm2)
00262 {
00263    REPORT
00264    int nr = gm->Nrows();
00265    MatrixRow mr(gm, LoadOnEntry+StoreOnExit+DirectPart);
00266    MatrixRow mr2(gm2, LoadOnEntry);
00267    while (nr--) { mr.RevSub(mr2); mr2.Next(); mr.Next(); }
00268 }
00269 
00270 static void SPDS(GeneralMatrix* gm, GeneralMatrix* gm1, GeneralMatrix* gm2)
00271 {
00272    REPORT
00273    int nr = gm->Nrows();
00274    MatrixRow mr1(gm1, LoadOnEntry); MatrixRow mr2(gm2, LoadOnEntry);
00275    MatrixRow mr(gm, StoreOnExit+DirectPart);
00276    while (nr--) { mr.Multiply(mr1,mr2); mr1.Next(); mr2.Next(); mr.Next(); }
00277 }
00278 
00279 static void SPDS(GeneralMatrix* gm, GeneralMatrix* gm2)
00280 // SP into first argument
00281 {
00282    REPORT
00283    int nr = gm->Nrows();
00284    MatrixRow mr(gm, StoreOnExit+LoadOnEntry+DirectPart);
00285    MatrixRow mr2(gm2, LoadOnEntry);
00286    while (nr--) { mr.Multiply(mr2); mr.Next(); mr2.Next(); }
00287 }
00288 
00289 static GeneralMatrix* GeneralMult1(GeneralMatrix* gm1, GeneralMatrix* gm2,
00290    MultipliedMatrix* mm, MatrixType mtx)
00291 {
00292    REPORT
00293    Tracer tr("GeneralMult1");
00294    int nr=gm1->Nrows(); int nc=gm2->Ncols();
00295    if (gm1->Ncols() !=gm2->Nrows())
00296       Throw(IncompatibleDimensionsException(*gm1, *gm2));
00297    GeneralMatrix* gmx = mtx.New(nr,nc,mm);
00298 
00299    MatrixCol mcx(gmx, StoreOnExit+DirectPart);
00300    MatrixCol mc2(gm2, LoadOnEntry);
00301    while (nc--)
00302    {
00303       MatrixRow mr1(gm1, LoadOnEntry, mcx.Skip());
00304       Real* el = mcx.Data();                         // pointer to an element
00305       int n = mcx.Storage();
00306       while (n--) { *(el++) = DotProd(mr1,mc2); mr1.Next(); }
00307       mc2.Next(); mcx.Next();
00308    }
00309    gmx->ReleaseAndDelete(); gm1->tDelete(); gm2->tDelete(); return gmx;
00310 }
00311 
00312 static GeneralMatrix* GeneralMult2(GeneralMatrix* gm1, GeneralMatrix* gm2,
00313    MultipliedMatrix* mm, MatrixType mtx)
00314 {
00315    // version that accesses by row only - not good for thin matrices
00316    // or column vectors in right hand term.
00317    REPORT
00318    Tracer tr("GeneralMult2");
00319    int nr=gm1->Nrows(); int nc=gm2->Ncols();
00320    if (gm1->Ncols() !=gm2->Nrows())
00321       Throw(IncompatibleDimensionsException(*gm1, *gm2));
00322    GeneralMatrix* gmx = mtx.New(nr,nc,mm);
00323 
00324    MatrixRow mrx(gmx, LoadOnEntry+StoreOnExit+DirectPart);
00325    MatrixRow mr1(gm1, LoadOnEntry);
00326    while (nr--)
00327    {
00328       MatrixRow mr2(gm2, LoadOnEntry, mr1.Skip());
00329       Real* el = mr1.Data();                         // pointer to an element
00330       int n = mr1.Storage();
00331       mrx.Zero();
00332       while (n--) { mrx.AddScaled(mr2, *el++); mr2.Next(); }
00333       mr1.Next(); mrx.Next();
00334    }
00335    gmx->ReleaseAndDelete(); gm1->tDelete(); gm2->tDelete(); return gmx;
00336 }
00337 
00338 static GeneralMatrix* mmMult(GeneralMatrix* gm1, GeneralMatrix* gm2)
00339 {
00340    // matrix multiplication for type Matrix only
00341    REPORT
00342    Tracer tr("MatrixMult");
00343 
00344    int nr=gm1->Nrows(); int ncr=gm1->Ncols(); int nc=gm2->Ncols();
00345    if (ncr != gm2->Nrows()) Throw(IncompatibleDimensionsException(*gm1,*gm2));
00346 
00347    Matrix* gm = new Matrix(nr,nc); MatrixErrorNoSpace(gm);
00348 
00349    Real* s1=gm1->Store(); Real* s2=gm2->Store(); Real* s=gm->Store();
00350 
00351    if (ncr)
00352    {
00353       while (nr--)
00354       {
00355          Real* s2x = s2; int j = ncr;
00356          Real* sx = s; Real f = *s1++; int k = nc;
00357          while (k--) *sx++ = f * *s2x++;
00358          while (--j)
00359             { sx = s; f = *s1++; k = nc; while (k--) *sx++ += f * *s2x++; }
00360          s = sx;
00361       }
00362    }
00363    else *gm = 0.0;
00364 
00365    gm->ReleaseAndDelete(); gm1->tDelete(); gm2->tDelete(); return gm;
00366 }
00367 
00368 static GeneralMatrix* GeneralMult(GeneralMatrix* gm1, GeneralMatrix* gm2,
00369    MultipliedMatrix* mm, MatrixType mtx)
00370 {
00371    if ( Rectangular(gm1->Type(), gm2->Type(), mtx))
00372    {
00373       REPORT
00374       return mmMult(gm1, gm2);
00375    }
00376    else
00377    {
00378       REPORT
00379       Compare(gm1->Type() * gm2->Type(),mtx);
00380       int nr = gm2->Nrows(); int nc = gm2->Ncols();
00381       if (nc <= 5 && nr > nc) { REPORT return GeneralMult1(gm1, gm2, mm, mtx); }
00382       else { REPORT return GeneralMult2(gm1, gm2, mm, mtx); }
00383    }
00384 }
00385 
00386 static GeneralMatrix* GeneralKP(GeneralMatrix* gm1, GeneralMatrix* gm2,
00387    KPMatrix* kp, MatrixType mtx)
00388 {
00389    REPORT
00390    Tracer tr("GeneralKP");
00391    int nr1 = gm1->Nrows(); int nc1 = gm1->Ncols();
00392    int nr2 = gm2->Nrows(); int nc2 = gm2->Ncols();
00393    Compare((gm1->Type()).KP(gm2->Type()),mtx);
00394    GeneralMatrix* gmx = mtx.New(nr1*nr2, nc1*nc2, kp);
00395    MatrixRow mrx(gmx, LoadOnEntry+StoreOnExit+DirectPart);
00396    MatrixRow mr1(gm1, LoadOnEntry);
00397    for (int i = 1; i <= nr1; ++i)
00398    {
00399       MatrixRow mr2(gm2, LoadOnEntry);
00400       for (int j = 1; j <= nr2; ++j)
00401          { mrx.KP(mr1,mr2); mr2.Next(); mrx.Next(); }
00402       mr1.Next();
00403    }
00404    gmx->ReleaseAndDelete(); gm1->tDelete(); gm2->tDelete(); return gmx;
00405 }
00406 
00407 static GeneralMatrix* GeneralSolv(GeneralMatrix* gm1, GeneralMatrix* gm2,
00408    BaseMatrix* sm, MatrixType mtx)
00409 {
00410    REPORT
00411    Tracer tr("GeneralSolv");
00412    Compare(gm1->Type().i() * gm2->Type(),mtx);
00413    int nr = gm1->Nrows();
00414    if (nr != gm1->Ncols()) Throw(NotSquareException(*gm1));
00415    int nc = gm2->Ncols();
00416    if (gm1->Ncols() != gm2->Nrows())
00417       Throw(IncompatibleDimensionsException(*gm1, *gm2));
00418    GeneralMatrix* gmx = mtx.New(nr,nc,sm); MatrixErrorNoSpace(gmx);
00419    Real* r = new Real [nr]; MatrixErrorNoSpace(r);
00420    MONITOR_REAL_NEW("Make   (GenSolv)",nr,r)
00421    GeneralMatrix* gms = gm1->MakeSolver();
00422    Try
00423    {
00424 
00425       MatrixColX mcx(gmx, r, StoreOnExit+DirectPart);   // copy to and from r
00426          // this must be inside Try so mcx is destroyed before gmx
00427       MatrixColX mc2(gm2, r, LoadOnEntry);
00428       while (nc--) { gms->Solver(mcx, mc2); mcx.Next(); mc2.Next(); }
00429    }
00430    CatchAll
00431    {
00432       if (gms) gms->tDelete();
00433       delete gmx;                   // <--------------------
00434       gm2->tDelete();
00435       MONITOR_REAL_DELETE("Delete (GenSolv)",nr,r)
00436                           // AT&T version 2.1 gives an internal error
00437       delete [] r;
00438       ReThrow;
00439    }
00440    gms->tDelete(); gmx->ReleaseAndDelete(); gm2->tDelete();
00441    MONITOR_REAL_DELETE("Delete (GenSolv)",nr,r)
00442                           // AT&T version 2.1 gives an internal error
00443    delete [] r;
00444    return gmx;
00445 }
00446 
00447 // version for inverses - gm2 is identity
00448 static GeneralMatrix* GeneralSolvI(GeneralMatrix* gm1, BaseMatrix* sm,
00449    MatrixType mtx)
00450 {
00451    REPORT
00452    Tracer tr("GeneralSolvI");
00453    Compare(gm1->Type().i(),mtx);
00454    int nr = gm1->Nrows();
00455    if (nr != gm1->Ncols()) Throw(NotSquareException(*gm1));
00456    int nc = nr;
00457    // DiagonalMatrix I(nr); I = 1;
00458    IdentityMatrix I(nr);
00459    GeneralMatrix* gmx = mtx.New(nr,nc,sm); MatrixErrorNoSpace(gmx);
00460    Real* r = new Real [nr]; MatrixErrorNoSpace(r);
00461    MONITOR_REAL_NEW("Make   (GenSolvI)",nr,r)
00462    GeneralMatrix* gms = gm1->MakeSolver();
00463    Try
00464    {
00465 
00466       MatrixColX mcx(gmx, r, StoreOnExit+DirectPart);   // copy to and from r
00467          // this must be inside Try so mcx is destroyed before gmx
00468       MatrixColX mc2(&I, r, LoadOnEntry);
00469       while (nc--) { gms->Solver(mcx, mc2); mcx.Next(); mc2.Next(); }
00470    }
00471    CatchAll
00472    {
00473       if (gms) gms->tDelete();
00474       delete gmx;
00475       MONITOR_REAL_DELETE("Delete (GenSolvI)",nr,r)
00476                           // AT&T version 2.1 gives an internal error
00477       delete [] r;
00478       ReThrow;
00479    }
00480    gms->tDelete(); gmx->ReleaseAndDelete();
00481    MONITOR_REAL_DELETE("Delete (GenSolvI)",nr,r)
00482                           // AT&T version 2.1 gives an internal error
00483    delete [] r;
00484    return gmx;
00485 }
00486 
00487 GeneralMatrix* InvertedMatrix::Evaluate(MatrixType mtx)
00488 {
00489    // Matrix Inversion - use solve routines
00490    Tracer tr("InvertedMatrix::Evaluate");
00491    REPORT
00492    gm=((BaseMatrix*&)bm)->Evaluate();
00493    return GeneralSolvI(gm,this,mtx);
00494 }
00495 
00496 //*************************** New versions ************************
00497 
00498 GeneralMatrix* AddedMatrix::Evaluate(MatrixType mtd)
00499 {
00500    REPORT
00501    Tracer tr("AddedMatrix::Evaluate");
00502    gm1=((BaseMatrix*&)bm1)->Evaluate(); gm2=((BaseMatrix*&)bm2)->Evaluate();
00503    int nr=gm1->Nrows(); int nc=gm1->Ncols();
00504    if (nr!=gm2->Nrows() || nc!=gm2->Ncols())
00505    {
00506       Try { Throw(IncompatibleDimensionsException(*gm1, *gm2)); }
00507       CatchAll
00508       {
00509          gm1->tDelete(); gm2->tDelete();
00510          ReThrow;
00511       }
00512    }
00513    MatrixType mt1 = gm1->Type(), mt2 = gm2->Type(); MatrixType mts = mt1 + mt2;
00514    if (!mtd) { REPORT mtd = mts; }
00515    else if (!(mtd.DataLossOK || mtd >= mts))
00516    {
00517       REPORT
00518       gm1->tDelete(); gm2->tDelete();
00519       Throw(ProgramException("Illegal Conversion", mts, mtd));
00520    }
00521    GeneralMatrix* gmx;
00522    bool c1 = (mtd == mt1), c2 = (mtd == mt2);
00523    if ( c1 && c2 && (gm1->SimpleAddOK(gm2) == 0) )
00524    {
00525       if (gm1->reuse()) { REPORT AddTo(gm1,gm2); gm2->tDelete(); gmx = gm1; }
00526       else if (gm2->reuse()) { REPORT AddTo(gm2,gm1); gmx = gm2; }
00527       else
00528       {
00529          REPORT
00530          // what if new throws an exception
00531          Try { gmx = mt1.New(nr,nc,this); }
00532          CatchAll
00533          {
00534             ReThrow;
00535          }
00536          gmx->ReleaseAndDelete(); Add(gmx,gm1,gm2);
00537       }
00538    }
00539    else
00540    {
00541       if (c1 && c2)
00542       {
00543          short SAO = gm1->SimpleAddOK(gm2);
00544          if (SAO & 1) { REPORT c1 = false; }
00545          if (SAO & 2) { REPORT c2 = false; }
00546       }
00547       if (c1 && gm1->reuse() )               // must have type test first
00548          { REPORT AddDS(gm1,gm2); gm2->tDelete(); gmx = gm1; }
00549       else if (c2 && gm2->reuse() )
00550          { REPORT AddDS(gm2,gm1); if (!c1) gm1->tDelete(); gmx = gm2; }
00551       else
00552       {
00553          REPORT
00554          Try { gmx = mtd.New(nr,nc,this); }
00555          CatchAll
00556          {
00557             if (!c1) gm1->tDelete(); if (!c2) gm2->tDelete();
00558             ReThrow;
00559          }
00560          AddDS(gmx,gm1,gm2);
00561          if (!c1) gm1->tDelete(); if (!c2) gm2->tDelete();
00562          gmx->ReleaseAndDelete();
00563       }
00564    }
00565    return gmx;
00566 }
00567 
00568 GeneralMatrix* SubtractedMatrix::Evaluate(MatrixType mtd)
00569 {
00570    REPORT
00571    Tracer tr("SubtractedMatrix::Evaluate");
00572    gm1=((BaseMatrix*&)bm1)->Evaluate(); gm2=((BaseMatrix*&)bm2)->Evaluate();
00573    int nr=gm1->Nrows(); int nc=gm1->Ncols();
00574    if (nr!=gm2->Nrows() || nc!=gm2->Ncols())
00575    {
00576       Try { Throw(IncompatibleDimensionsException(*gm1, *gm2)); }
00577       CatchAll
00578       {
00579          gm1->tDelete(); gm2->tDelete();
00580          ReThrow;
00581       }
00582    }
00583    MatrixType mt1 = gm1->Type(), mt2 = gm2->Type(); MatrixType mts = mt1 + mt2;
00584    if (!mtd) { REPORT mtd = mts; }
00585    else if (!(mtd.DataLossOK || mtd >= mts))
00586    {
00587       gm1->tDelete(); gm2->tDelete();
00588       Throw(ProgramException("Illegal Conversion", mts, mtd));
00589    }
00590    GeneralMatrix* gmx;
00591    bool c1 = (mtd == mt1), c2 = (mtd == mt2);
00592    if ( c1 && c2 && (gm1->SimpleAddOK(gm2) == 0) )
00593    {
00594       if (gm1->reuse())
00595          { REPORT SubtractFrom(gm1,gm2); gm2->tDelete(); gmx = gm1; }
00596       else if (gm2->reuse()) { REPORT ReverseSubtract(gm2,gm1); gmx = gm2; }
00597       else
00598       {
00599          REPORT
00600          Try { gmx = mt1.New(nr,nc,this); }
00601          CatchAll
00602          {
00603             ReThrow;
00604          }
00605          gmx->ReleaseAndDelete(); Subtract(gmx,gm1,gm2);
00606       }
00607    }
00608    else
00609    {
00610       if (c1 && c2)
00611       {
00612          short SAO = gm1->SimpleAddOK(gm2);
00613          if (SAO & 1) { REPORT c1 = false; }
00614          if (SAO & 2) { REPORT c2 = false; }
00615       }
00616       if (c1 && gm1->reuse() )               // must have type test first
00617          { REPORT SubtractDS(gm1,gm2); gm2->tDelete(); gmx = gm1; }
00618       else if (c2 && gm2->reuse() )
00619       {
00620          REPORT ReverseSubtractDS(gm2,gm1);
00621          if (!c1) gm1->tDelete(); gmx = gm2;
00622       }
00623       else
00624       {
00625          REPORT
00626          // what if New throws and exception
00627          Try { gmx = mtd.New(nr,nc,this); }
00628          CatchAll
00629          {
00630             if (!c1) gm1->tDelete(); if (!c2) gm2->tDelete();
00631             ReThrow;
00632          }
00633          SubtractDS(gmx,gm1,gm2);
00634          if (!c1) gm1->tDelete(); if (!c2) gm2->tDelete();
00635          gmx->ReleaseAndDelete();
00636       }
00637    }
00638    return gmx;
00639 }
00640 
00641 GeneralMatrix* SPMatrix::Evaluate(MatrixType mtd)
00642 {
00643    REPORT
00644    Tracer tr("SPMatrix::Evaluate");
00645    gm1=((BaseMatrix*&)bm1)->Evaluate(); gm2=((BaseMatrix*&)bm2)->Evaluate();
00646    int nr=gm1->Nrows(); int nc=gm1->Ncols();
00647    if (nr!=gm2->Nrows() || nc!=gm2->Ncols())
00648    {
00649       Try { Throw(IncompatibleDimensionsException(*gm1, *gm2)); }
00650       CatchAll
00651       {
00652          gm1->tDelete(); gm2->tDelete();
00653          ReThrow;
00654       }
00655    }
00656    MatrixType mt1 = gm1->Type(), mt2 = gm2->Type();
00657    MatrixType mts = mt1.SP(mt2);
00658    if (!mtd) { REPORT mtd = mts; }
00659    else if (!(mtd.DataLossOK || mtd >= mts))
00660    {
00661       gm1->tDelete(); gm2->tDelete();
00662       Throw(ProgramException("Illegal Conversion", mts, mtd));
00663    }
00664    GeneralMatrix* gmx;
00665    bool c1 = (mtd == mt1), c2 = (mtd == mt2);
00666    if ( c1 && c2 && (gm1->SimpleAddOK(gm2) == 0) )
00667    {
00668       if (gm1->reuse()) { REPORT SP(gm1,gm2); gm2->tDelete(); gmx = gm1; }
00669       else if (gm2->reuse()) { REPORT SP(gm2,gm1); gmx = gm2; }
00670       else
00671       {
00672          REPORT
00673          Try { gmx = mt1.New(nr,nc,this); }
00674          CatchAll
00675          {
00676             ReThrow;
00677          }
00678          gmx->ReleaseAndDelete(); SP(gmx,gm1,gm2);
00679       }
00680    }
00681    else
00682    {
00683       if (c1 && c2)
00684       {
00685          short SAO = gm1->SimpleAddOK(gm2);
00686          if (SAO & 1) { REPORT c2 = false; }    // c1 and c2 swapped
00687          if (SAO & 2) { REPORT c1 = false; }
00688       }
00689       if (c1 && gm1->reuse() )               // must have type test first
00690          { REPORT SPDS(gm1,gm2); gm2->tDelete(); gmx = gm1; }
00691       else if (c2 && gm2->reuse() )
00692          { REPORT SPDS(gm2,gm1); if (!c1) gm1->tDelete(); gmx = gm2; }
00693       else
00694       {
00695          REPORT
00696          // what if New throws and exception
00697          Try { gmx = mtd.New(nr,nc,this); }
00698          CatchAll
00699          {
00700             if (!c1) gm1->tDelete(); if (!c2) gm2->tDelete();
00701             ReThrow;
00702          }
00703          SPDS(gmx,gm1,gm2);
00704          if (!c1) gm1->tDelete(); if (!c2) gm2->tDelete();
00705          gmx->ReleaseAndDelete();
00706       }
00707    }
00708    return gmx;
00709 }
00710 
00711 
00712 
00713 //*************************** norm functions ****************************/
00714 
00715 Real BaseMatrix::Norm1() const
00716 {
00717    // maximum of sum of absolute values of a column
00718    REPORT
00719    GeneralMatrix* gm = ((BaseMatrix&)*this).Evaluate();
00720    int nc = gm->Ncols(); Real value = 0.0;
00721    MatrixCol mc(gm, LoadOnEntry);
00722    while (nc--)
00723       { Real v = mc.SumAbsoluteValue(); if (value < v) value = v; mc.Next(); }
00724    gm->tDelete(); return value;
00725 }
00726 
00727 Real BaseMatrix::NormInfinity() const
00728 {
00729    // maximum of sum of absolute values of a row
00730    REPORT
00731    GeneralMatrix* gm = ((BaseMatrix&)*this).Evaluate();
00732    int nr = gm->Nrows(); Real value = 0.0;
00733    MatrixRow mr(gm, LoadOnEntry);
00734    while (nr--)
00735       { Real v = mr.SumAbsoluteValue(); if (value < v) value = v; mr.Next(); }
00736    gm->tDelete(); return value;
00737 }
00738 
00739 //********************** Concatenation and stacking *************************/
00740 
00741 GeneralMatrix* ConcatenatedMatrix::Evaluate(MatrixType mtx)
00742 {
00743    REPORT
00744    Tracer tr("Concatenate");
00745       gm2 = ((BaseMatrix*&)bm2)->Evaluate();
00746       gm1 = ((BaseMatrix*&)bm1)->Evaluate();
00747       Compare(gm1->Type() | gm2->Type(),mtx);
00748       int nr=gm1->Nrows(); int nc = gm1->Ncols() + gm2->Ncols();
00749       if (nr != gm2->Nrows())
00750          Throw(IncompatibleDimensionsException(*gm1, *gm2));
00751