My title

00001 /* ****************************************************************** **
00002 **    OpenSees - Open System for Earthquake Engineering Simulation    **
00003 **          Pacific Earthquake Engineering Research Center            **
00004 **                                                                    **
00005 **                                                                    **
00006 ** (C) Copyright 1999, The Regents of the University of California    **
00007 ** All Rights Reserved.                                               **
00008 **                                                                    **
00009 ** Commercial use of this program without express permission of the   **
00010 ** University of California, Berkeley, is strictly prohibited.  See   **
00011 ** file 'COPYRIGHT'  in main directory for information on usage and   **
00012 ** redistribution,  and for a DISCLAIMER OF ALL WARRANTIES.           **
00013 **                                                                    **
00014 ** Developed by:                                                      **
00015 **   Frank McKenna (fmckenna@ce.berkeley.edu)                         **
00016 **   Gregory L. Fenves (fenves@ce.berkeley.edu)                       **
00017 **   Filip C. Filippou (filippou@ce.berkeley.edu)                     **
00018 **                                                                    **
00019 ** ****************************************************************** */
00020                                                                         
00021 // $Revision: 1.4 $
00022 // $Date: 2001/02/17 04:03:10 $
00023 // $Source: /usr/local/cvs/OpenSees/SRC/matrix/Matrix.cpp,v $
00024                                                                         
00025                                                                         
00026 // File: ~/matrix/Matrix.h
00027 //
00028 // Written: fmk 
00029 // Created: 11/96
00030 // Revision: A
00031 //
00032 // Description: This file contains the class implementation for Matrix.
00033 //
00034 // What: "@(#) Matrix.h, revA"
00035 
00036 #include "Matrix.h"
00037 #include "Vector.h"
00038 #include "ID.h"
00039 
00040 #include <stdlib.h>
00041 
00042 #define MATRIX_WORK_AREA 400
00043 #define INT_WORK_AREA 20
00044 
00045 double Matrix::MATRIX_NOT_VALID_ENTRY =0.0;
00046 double *Matrix::matrixWork = new double [MATRIX_WORK_AREA];
00047 int *Matrix::intWork = new int [INT_WORK_AREA];
00048 //double *Matrix::matrixWork = (double *)malloc(400*sizeof(double));
00049 
00050 //
00051 // CONSTRUCTORS
00052 //
00053 
00054 Matrix::Matrix()
00055 :numRows(0), numCols(0),dataSize(0),data(0),fromFree(0)
00056 {
00057     // does nothing
00058 }
00059 
00060 
00061 Matrix::Matrix(int nRows,int nCols)
00062 :numRows(nRows),numCols(nCols),fromFree(0)
00063 {
00064 #ifdef _G3DEBUG
00065     if (nRows < 0) {
00066       cerr << "WARNING: Matrix::Matrix(int,int): tried to init matrix ";
00067       cerr << "with num rows: " << nRows << " <0\n";
00068       numRows = 0; numCols =0; dataSize =0; data = 0;
00069     }
00070     if (nCols < 0) {
00071       cerr << "WARNING: Matrix::Matrix(int,int): tried to init matrix";
00072       cerr << "with num cols: " << nCols << " <0\n";
00073       numRows = 0; numCols =0; dataSize =0; data = 0;
00074     }
00075 #endif
00076     dataSize = numRows * numCols;
00077     data = 0;
00078 
00079     if (dataSize > 0) {
00080       data = new double[dataSize];
00081       //data = (double *)malloc(dataSize*sizeof(double));
00082       if (data == 0) {
00083  cerr << "WARNING:Matrix::Matrix(int,int): Ran out of memory on init ";
00084  cerr << "of size " << dataSize << "\n";
00085  numRows = 0; numCols =0; dataSize =0;
00086       } else {
00087  // zero the data
00088  double *dataPtr = data;
00089  for (int i=0; i<dataSize; i++)
00090    *dataPtr++ = 0.0;
00091       }
00092     }
00093 }
00094 
00095 Matrix::Matrix(double *theData, int row, int col) 
00096 :numRows(row),numCols(col),dataSize(row*col),data(theData),fromFree(1)
00097 {
00098 #ifdef _G3DEBUG
00099     if (row < 0) {
00100       cerr << "WARNING: Matrix::Matrix(int,int): tried to init matrix with numRows: ";
00101       cerr << row << " <0\n";
00102       numRows = 0; numCols =0; dataSize =0; data = 0;
00103     }
00104     if (col < 0) {
00105       cerr << "WARNING: Matrix::Matrix(int,int): tried to init matrix with numCols: ";
00106       cerr << col << " <0\n";
00107       numRows = 0; numCols =0; dataSize =0; data = 0;
00108     }    
00109 #endif
00110 
00111     // does nothing
00112 }
00113 
00114 Matrix::Matrix(const Matrix &other)
00115 :fromFree(0)
00116 {
00117     numRows = other.numRows;
00118     numCols = other.numCols;
00119     dataSize = other.dataSize;
00120     data = 0;
00121     
00122     if (dataSize != 0) {
00123       data = new double[dataSize];
00124       // data = (double *)malloc(dataSize*sizeof(double));
00125       if (data == 0) {
00126  cerr << "WARNING:Matrix::Matrix(Matrix &): ";
00127  cerr << "Ran out of memory on init of size " << dataSize << "\n"; 
00128  numRows = 0; numCols =0; dataSize = 0;
00129       } else {
00130  // copy the data
00131  double *dataPtr = data;
00132  double *otherDataPtr = other.data;
00133  for (int i=0; i<dataSize; i++)
00134    *dataPtr++ = *otherDataPtr++;
00135       }
00136     }
00137 }
00138 
00139 
00140 //
00141 // DESTRUCTOR
00142 //
00143 
00144 Matrix::~Matrix()
00145 {
00146   if (data != 0) 
00147       if (fromFree == 0)
00148    delete [] data; 
00149   //  if (data != 0) free((void *) data);
00150 }
00151     
00152 
00153 //
00154 // METHODS - Zero, Assemble, Solve
00155 //
00156 
00157 void
00158 Matrix::Zero(void)
00159 {
00160   double *dataPtr = data;
00161   for (int i=0; i<dataSize; i++)
00162     *dataPtr++ = 0;
00163 }
00164 
00165 
00166 int
00167 Matrix::resize(int rows, int cols) {
00168 
00169   int newSize = rows*cols;
00170 
00171   if (newSize <= 0) {
00172     g3ErrorHandler->warning("Matrix::resize) - rows %d or cols %d specified <= 0\n",
00173        rows, cols); 
00174     return -1;
00175   }
00176 
00177   else if (newSize > dataSize) {
00178 
00179     // free the old space
00180     if (data != 0) 
00181       if (fromFree == 0)
00182  delete [] data; 
00183     //  if (data != 0) free((void *) data);
00184 
00185     fromFree = 0;
00186     // create new space
00187     data = new double[newSize];
00188     // data = (double *)malloc(dataSize*sizeof(double));
00189     if (data == 0) {
00190       g3ErrorHandler->warning("Matrix::resize(%d %d) - out of memory\n",
00191          rows, cols); 
00192       numRows = 0; numCols =0; dataSize = 0;
00193       return -2;
00194     }
00195     dataSize = newSize;
00196     numRows = rows;
00197     numCols = cols;
00198   }
00199 
00200   // just reset the cols and rows - save two memory calls at expense of holding 
00201   // onto extra memory
00202   else {
00203     numRows = rows;
00204     numCols = cols;
00205   }
00206 
00207   return 0;
00208 }
00209 
00210 
00211 
00212 
00213 
00214 int
00215 Matrix::Assemble(const Matrix &V, const ID &rows, const ID &cols, double fact) 
00216 {
00217   int pos_Rows, pos_Cols;
00218   int res = 0;
00219 
00220   for (int i=0; i<cols.Size(); i++) {
00221     pos_Cols = cols(i);
00222     for (int j=0; j<rows.Size(); j++) {
00223       pos_Rows = rows(j);
00224       
00225       if ((pos_Cols >= 0) && (pos_Rows >= 0) && (pos_Rows < numRows) &&
00226    (pos_Cols < numCols) && (i < V.numCols) && (j < V.numRows))
00227  (*this)(pos_Rows,pos_Cols) += V(j,i)*fact;
00228       else {
00229  cerr << "WARNING: Matrix::Assemble(const Matrix &V, const ID &l): ";
00230  cerr << " - position (" << pos_Rows << "," << pos_Cols << ") outside bounds \n";
00231  res = -1;
00232       }
00233     }
00234   }
00235 
00236   return res;
00237 }
00238 
00239 #ifdef _WIN32
00240 extern "C" int _stdcall DGESV(int *N, int *NRHS, double *A, int *LDA, 
00241          int *iPiv, double *B, int *LDB, int *INFO);
00242         
00243 extern "C" int _stdcall DGETRS(char *TRANS, unsigned int sizeT,
00244           int *N, int *NRHS, double *A, int *LDA, 
00245           int *iPiv, double *B, int *LDB, int *INFO);
00246 #else
00247 extern "C" int dgesv_(int *N, int *NRHS, double *A, int *LDA, int *iPiv, 
00248         double *B, int *LDB, int *INFO);
00249 
00250 extern "C" int dgetrs_(char *TRANS, int *N, int *NRHS, double *A, int *LDA, 
00251          int *iPiv, double *B, int *LDB, int *INFO);         
00252 #endif
00253 
00254 int
00255 Matrix::Solve(const Vector &b, Vector &x) const
00256 {
00257 
00258     int n = numRows;
00259 
00260 #ifdef _G3DEBUG    
00261     if (numRows != numCols) {
00262       g3ErrorHandler->warning("Matrix::Solve(b,x) - %s [%d %d] is not square\n",
00263          "the matrix of dimensions",
00264          numRows, numCols);
00265       return -1;
00266     }
00267 
00268     if (n != x.Size()) {
00269       g3ErrorHandler->warning("Matrix::Solve(b,x) - dimension of x, %d, %s %d\n",
00270          numRows, "is not same as matrix", x.Size());
00271       return -2;
00272     }
00273 
00274     if (n != b.Size()) {
00275       g3ErrorHandler->warning("Matrix::Solve(b,x) - dimension of b, %d, %s %d\n",
00276          numRows, "is not same as matrix", b.Size());
00277       return -2;
00278     }
00279 #endif
00280     
00281     // check work area can hold all the data
00282     if (dataSize > MATRIX_WORK_AREA) {
00283        g3ErrorHandler->warning("Matrix::Solve(b,x) - matrix dimension [%d %d] %s %d\n",
00284           numRows, numCols, "larger than work area",
00285           MATRIX_WORK_AREA);
00286        return -3;
00287     }
00288 
00289     // copy the data
00290     int i;
00291     for (i=0; i<dataSize; i++)
00292       matrixWork[i] = data[i];
00293 
00294     // set x equal to b
00295     x = b;
00296 
00297     int nrhs = 1;
00298     int ldA = n;
00299     int ldB = n;
00300     int info;
00301     double *Aptr = matrixWork;
00302     double *Xptr = x.theData;
00303     int *iPIV = intWork;
00304     
00305 
00306 #ifdef _WIN32
00307     DGESV(&n,&nrhs,Aptr,&ldA,iPIV,Xptr,&ldB,&info);
00308 #else
00309     dgesv_(&n,&nrhs,Aptr,&ldA,iPIV,Xptr,&ldB,&info);
00310 #endif
00311 
00312 
00313     return 0;
00314 }
00315 
00316 
00317 int
00318 Matrix::Solve(const Matrix &b, Matrix &x) const
00319 {
00320 
00321     int n = numRows;
00322     int nrhs = x.numCols;
00323 
00324 #ifdef _G3DEBUG    
00325     if (numRows != numCols) {
00326       g3ErrorHandler->warning("Matrix::Solve(B,X) - the matrix of dimensions [%d %d] is not square\n",
00327          numRows, numCols);
00328       return -1;
00329     }
00330 
00331     if (n != x.numRows) {
00332       g3ErrorHandler->warning("Matrix::Solve(B,X) - #rows of X, %d, is not same as matrix %d\n",
00333          numRows, x.numRows);
00334       return -2;
00335     }
00336 
00337     if (n != b.numRows) {
00338       g3ErrorHandler->warning("Matrix::Solve(B,X) - #rows of B, %d, is not same as matrix %d\n",
00339          numRows, b.numRows);
00340       return -2;
00341     }
00342 
00343     if (x.numCols != b.numCols) {
00344       g3ErrorHandler->warning("Matrix::Solve(B,X) - #cols of B, %d, is not same as that of X, %d\n",
00345          b.numCols, x.numCols);
00346       return -3;
00347     }
00348 #endif
00349     
00350     // check work area can hold all the data
00351     if (dataSize > MATRIX_WORK_AREA || n > INT_WORK_AREA) {
00352       g3ErrorHandler->warning("Matrix::Solve(b,x) - matrix dimension [%d %d] larger than work area %d\n",
00353          numRows, numCols, MATRIX_WORK_AREA);
00354       return -3;
00355     }
00356 
00357     x = b;
00358 
00359     // copy the data
00360     int i;
00361     for (i=0; i<dataSize; i++)
00362       matrixWork[i] = data[i];
00363 
00364 
00365     int ldA = n;
00366     int ldB = n;
00367     int info;
00368     double *Aptr = matrixWork;
00369     double *Xptr = x.data;
00370     
00371     int *iPIV = intWork;
00372     
00373 
00374 #ifdef _WIN32
00375     DGESV(&n,&nrhs,Aptr,&ldA,iPIV,Xptr,&ldB,&info);
00376 #else
00377     dgesv_(&n,&nrhs,Aptr,&ldA,iPIV,Xptr,&ldB,&info);
00378 #endif
00379 
00380     return info;
00381 }
00382     
00383       
00384 
00385 int
00386 Matrix::addMatrix(double factThis, const Matrix &other, double factOther)
00387 {
00388     if (factThis == 1.0 && factOther == 0.0)
00389       return 0;
00390 
00391 #ifdef _G3DEBUG
00392     if ((other.numRows != numRows) || (other.numCols != numCols)) {
00393       g3ErrorHandler->warning("Matrix::addMatrix(): incompatable matrices, this[%d %d] other[%d %d]\n",
00394          numRows, numCols, other.numRows, other.numCols);
00395       return -1;
00396     }
00397 #endif
00398 
00399     if (factThis == 1.0) {
00400 
00401       // want: this += other * factOther
00402       if (factOther == 1.0) {
00403  double *dataPtr = data;
00404  double *otherDataPtr = other.data;      
00405  for (int i=0; i<dataSize; i++)
00406    *dataPtr++ += *otherDataPtr++;
00407       } else {
00408  double *dataPtr = data;
00409  double *otherDataPtr = other.data;      
00410  for (int i=0; i<dataSize; i++)
00411    *dataPtr++ += *otherDataPtr++ * factOther;
00412       }
00413     } 
00414 
00415     else if (factThis == 0.0) {
00416 
00417       // want: this = other * factOther
00418       if (factOther == 1.0) {
00419  double *dataPtr = data;
00420  double *otherDataPtr = other.data;      
00421  for (int i=0; i<dataSize; i++)
00422    *dataPtr++ = *otherDataPtr++;
00423       } else {
00424  double *dataPtr = data;
00425  double *otherDataPtr = other.data;      
00426  for (int i=0; i<dataSize; i++)
00427    *dataPtr++ = *otherDataPtr++ * factOther;
00428       }
00429     } 
00430 
00431     else {
00432 
00433       // want: this = this * thisFact + other * factOther
00434       if (factOther == 1.0) {
00435  double *dataPtr = data;
00436  double *otherDataPtr = other.data;      
00437  for (int i=0; i<dataSize; i++) {
00438    double value = *dataPtr * factThis + *otherDataPtr++;
00439    *dataPtr++ = value;
00440  }
00441       } else {
00442  double *dataPtr = data;
00443  double *otherDataPtr = other.data;      
00444  for (int i=0; i<dataSize; i++) {
00445    double value = *dataPtr * factThis + *otherDataPtr++ * factOther;
00446    *dataPtr++ = value;
00447  }
00448       }
00449     } 
00450 
00451     // successfull
00452     return 0;
00453 }
00454 
00455 
00456 
00457 int
00458 Matrix::addMatrixProduct(double thisFact, 
00459     const Matrix &B, 
00460     const Matrix &C, 
00461     double otherFact)
00462 {
00463     if (thisFact == 1.0 && otherFact == 0.0)
00464       return 0;
00465 #ifdef _G3DEBUG
00466     if ((B.numRows != numRows) || (C.numCols != numCols) || (B.numCols != C.numRows)) {
00467       g3ErrorHandler->warning("Matrix::addMatrixProduct():%s[%d %d] A[%d &d] B{%d %d]\n", 
00468          "incompatable matrices, this",
00469          numRows, numCols, B.numRows, B.numCols, C.numRows, C.numCols);
00470       return -1;
00471     }
00472 #endif
00473     // NOTE: looping as per blas3 dgemm_: j,k,i
00474     if (thisFact == 1.0) {
00475 
00476       if (otherFact == 1.0) {
00477  // want: this += B * C 
00478  int numColB = B.numCols;
00479  double *ckjPtr  = &(C.data)[0];
00480  for (int j=0; j<numCols; j++) {
00481    double *aijPtrA = &data[j*numRows];
00482    for (int k=0; k<numColB; k++) {
00483      double tmp = *ckjPtr++;
00484      double *aijPtr = aijPtrA;
00485      double *bikPtr = &(B.data)[k*numRows];
00486      for (int i=0; i<numRows; i++)
00487        *aijPtr++ += *bikPtr++ * tmp;
00488    }
00489  }
00490       } else {
00491  // want: this += B * C * otherFact
00492  int numColB = B.numCols;
00493  double *ckjPtr  = &(C.data)[0];
00494  for (int j=0; j<numCols; j++) {
00495    double *aijPtrA = &data[j*numRows];
00496    for (int k=0; k<numColB; k++) {
00497      double tmp = *ckjPtr++ * otherFact;
00498      double *aijPtr = aijPtrA;
00499      double *bikPtr = &(B.data)[k*numRows];
00500      for (int i=0; i<numRows; i++)
00501        *aijPtr++ += *bikPtr++ * tmp;
00502    }
00503  }
00504       }
00505     }
00506 
00507     else if (thisFact == 0.0) {
00508 
00509       if (otherFact == 1.0) {
00510  // want: this = B * C  otherFact
00511  double *dataPtr = data;
00512  for (int i=0; i<dataSize; i++)
00513    *dataPtr++ = 0.0;
00514  int numColB = B.numCols;
00515  double *ckjPtr  = &(C.data)[0];
00516  for (int j=0; j<numCols; j++) {
00517    double *aijPtrA = &data[j*numRows];
00518    for (int k=0; k<numColB; k++) {
00519      double tmp = *ckjPtr++;
00520      double *aijPtr = aijPtrA;
00521      double *bikPtr = &(B.data)[k*numRows];
00522      for (int i=0; i<numRows; i++)
00523        *aijPtr++ += *bikPtr++ * tmp;
00524    }
00525  }
00526       } else {
00527  // want: this = B * C  otherFact
00528  double *dataPtr = data;
00529  for (int i=0; i<dataSize; i++)
00530    *dataPtr++ = 0.0;
00531  int numColB = B.numCols;
00532  double *ckjPtr  = &(C.data)[0];
00533  for (int j=0; j<numCols; j++) {
00534    double *aijPtrA = &data[j*numRows];
00535    for (int k=0; k<numColB; k++) {
00536      double tmp = *ckjPtr++ * otherFact;
00537      double *aijPtr = aijPtrA;
00538      double *bikPtr = &(B.data)[k*numRows];
00539      for (int i=0; i<numRows; i++)
00540        *aijPtr++ += *bikPtr++ * tmp;
00541    }
00542  }
00543       }
00544     } 
00545 
00546     else {
00547       // want: this = thisFact * this + B * C  otherFact
00548       double *dataPtr = data;
00549       for (int i=0; i<dataSize; i++)
00550    *dataPtr++ *= thisFact;
00551       int numColB = B.numCols;
00552       double *ckjPtr  = &(C.data)[0];
00553       for (int j=0; j<numCols; j++) {
00554  double *aijPtrA = &data[j*numRows];
00555  for (int k=0; k<numColB; k++) {
00556    double tmp = *ckjPtr++ * otherFact;
00557    double *aijPtr = aijPtrA;
00558    double *bikPtr = &(B.data)[k*numRows];
00559    for (int i=0; i<numRows; i++)
00560      *aijPtr++ += *bikPtr++ * tmp;
00561  }
00562       }
00563     } 
00564 
00565     return 0;
00566 }
00567 
00568 
00569 // to perform this += T' * B * T
00570 int
00571 Matrix::addMatrixTripleProduct(double thisFact, 
00572           const Matrix &T, 
00573           const Matrix &B, 
00574           double otherFact)
00575 {
00576     if (thisFact == 1.0 && otherFact == 0.0)
00577       return 0;
00578 #ifdef _G3DEBUG
00579     if ((numCols != numRows) || (B.numCols != B.numRows) || (T.numCols != numRows) ||
00580  (T.numRows != B.numCols)) {
00581       g3ErrorHandler->warning("Matrix::addMatrixTripleProduct():%s[%d %d] T[%d %d] B[%d %d]\n", 
00582          "incompatable matrices, this",
00583          numRows, numCols, T.numRows, T.numCols, B.numRows, B.numCols);
00584       return -1;
00585     }
00586 #endif
00587 
00588     // cheack work area can hold the temporary matrix
00589     int dimB = B.numCols;
00590     int sizeWork = dimB * numCols;
00591 
00592     if (sizeWork > MATRIX_WORK_AREA) {
00593       this->addMatrix(thisFact, T^B*T, otherFact);
00594       return 0;
00595     }
00596 
00597     // zero out the work area
00598     double *matrixWorkPtr = matrixWork;
00599     for (int l=0; l<sizeWork; l++)
00600       *matrixWorkPtr++ = 0.0;
00601 
00602     // now form B * T * fact store in matrixWork == A area
00603     // NOTE: looping as per blas3 dgemm_: j,k,i
00604 
00605     double *tkjPtr  = &(T.data)[0];
00606     for (int j=0; j<numCols; j++) {
00607       double *aijPtrA = &matrixWork[j*dimB];
00608       for (int k=0; k<dimB; k++) {
00609  double tmp = *tkjPtr++ * otherFact;
00610  double *aijPtr = aijPtrA;
00611  double *bikPtr = &(B.data)[k*dimB];
00612  for (int i=0; i<dimB; i++) 
00613    *aijPtr++ += *bikPtr++ * tmp;
00614       }
00615     }
00616 
00617     // now form T' * matrixWork
00618     // NOTE: looping as per blas3 dgemm_: j,i,k
00619     if (thisFact == 1.0) {
00620       double *dataPtr = &data[0];
00621       for (int j=0; j< numCols; j++) {
00622  double *workkjPtrA = &matrixWork[j*dimB];
00623  for (int i=0; i<numRows; i++) {
00624    double *ckiPtr = &(T.data)[i*dimB];
00625    double *workkjPtr = workkjPtrA;
00626    double aij = 0.0;
00627    for (int k=0; k< dimB; k++)
00628      aij += *ckiPtr++ * *workkjPtr++;
00629    *dataPtr++ += aij;
00630  }
00631       }
00632     } else if (thisFact == 0.0) {
00633       double *dataPtr = &data[0];
00634       for (int j=0; j< numCols; j++) {
00635  double *workkjPtrA = &matrixWork[j*dimB];
00636  for (int i=0; i<numRows; i++) {
00637    double *ckiPtr = &(T.data)[i*dimB];
00638    double *workkjPtr = workkjPtrA;
00639    double aij = 0.0;
00640    for (int k=0; k< dimB; k++)
00641      aij += *ckiPtr++ * *workkjPtr++;
00642    *dataPtr++ = aij;
00643  }
00644       }
00645 
00646     } else {
00647       double *dataPtr = &data[0];
00648       for (int j=0; j< numCols; j++) {
00649  double *workkjPtrA = &matrixWork[j*dimB];
00650  for (int i=0; i<numRows; i++) {
00651    double *ckiPtr = &(T.data)[i*dimB];
00652    double *workkjPtr = workkjPtrA;
00653    double aij = 0.0;
00654    for (int k=0; k< dimB; k++)
00655      aij += *ckiPtr++ * *workkjPtr++;
00656    double value = *dataPtr * thisFact + aij;
00657    *dataPtr++ = value;
00658  }
00659       }
00660     }
00661 
00662     return 0;
00663 }
00664 
00665 
00666 
00667 //
00668 // OVERLOADED OPERATOR () to CONSTRUCT A NEW MATRIX
00669 //
00670 
00671 Matrix
00672 Matrix::operator()(const ID &rows, const ID & cols) const
00673 {
00674     int nRows, nCols;
00675     nRows = rows.Size();
00676     nCols = cols.Size();
00677     Matrix result(nRows,nCols);
00678     double *dataPtr = result.data;
00679     for (int i=0; i<nCols; i++)
00680  for (int j=0; j<nRows; j++)
00681      *dataPtr++ = (*this)(rows(j),cols(i));
00682 
00683     return result;
00684 }
00685   
00686 // Matrix &operator=(const Matrix  &V):
00687 //      the assignment operator, This is assigned to be a copy of V. if sizes
00688 //      are not compatable this.data [] is deleted. The data pointers will not
00689 //      point to the same area in mem after the assignment.
00690 //
00691 
00692 
00693 
00694 Matrix &
00695 Matrix::operator=(const Matrix &other)
00696 {
00697   // first check we are not trying other = other
00698   if (this == &other) 
00699     return *this;
00700 
00701 /*
00702 #ifdef _G3DEBUG    
00703   if ((numCols != other.numCols) || (numRows != other.numRows)) {
00704     g3ErrorHandler->warning("Matrix::operator=() - matrix dimensions do not match: [%d %d] != [%d %d]\n",
00705        numRows, numCols, other.numRows, other.numCols);
00706     return *this;
00707   }
00708 #endif
00709 */
00710   if ((numCols != other.numCols) || (numRows != other.numRows)) {
00711 #ifdef _G3DEBUG    
00712       g3ErrorHandler->warning("Matrix::operator=() - matrix dimensions do not match: [%d %d] != [%d %d]\n",
00713          numRows, numCols, other.numRows, other.numCols);
00714 #endif
00715       if (this->data != 0)
00716    delete [] this->data;
00717       
00718       int theSize = other.numCols*other.numRows;
00719       
00720       data = new double[theSize];
00721       
00722       this->dataSize = theSize;
00723       this->numCols = other.numCols;
00724       this->numRows = other.numRows;
00725   }
00726 
00727 
00728   // now copy the data
00729   double *dataPtr = data;
00730   double *otherDataPtr = other.data;      
00731   for (int i=0; i<dataSize; i++)
00732       *dataPtr++ = *otherDataPtr++;
00733   
00734   return *this;
00735 }
00736 
00737 
00738 
00739 
00740 
00741 Matrix &
00742 Matrix::operator=(const Tensor &V)
00743 {
00744   int rank = V.rank();
00745   if (rank != 4) {
00746       g3ErrorHandler->warning("Matrix::operator=() - tensor must be of rank 4\n");
00747       return *this;
00748   }
00749   int dim = V.dim(1);
00750   if (dim != V.dim(2) != V.dim(3) != V.dim(4)) {
00751       g3ErrorHandler->warning("Matrix::operator=() - tensor must have square dimensions\n");
00752       return *this;
00753   }
00754 
00755   if (dim != 2 || dim != 3 || dim != 1) {
00756       g3ErrorHandler->warning("Matrix::operator=() - tensor must be of dimension 2 or 3\n");
00757       return *this;
00758   }      
00759 
00760   if (dim == 1) {
00761       if ((numCols != 1) || (numRows != 1)) {      
00762    g3ErrorHandler->warning("Matrix::operator=() - matrix must be %s\n",
00763       "1x1 for tensor of dimension 3\n");
00764    return *this;
00765       }         
00766       (*this)(0,0) = V.cval(1,1,1,1);
00767       
00768   } else if (dim == 2) {
00769       if ((numCols != 3) || (numRows != 3)) {      
00770    g3ErrorHandler->warning("Matrix::operator=() - matrix must be %s\n",
00771       "1x1 for tensor of dimension 3\n");      
00772    return *this;
00773       }
00774       (*this)(0,0) = V.cval(1,1,1,1);
00775       (*this)(0,1) = V.cval(1,1,2,2);
00776       (*this)(0,2) = V.cval(1,1,1,2);      
00777 
00778       (*this)(1,0) = V.cval(2,2,1,1);
00779       (*this)(1,1) = V.cval(2,2,2,2);
00780       (*this)(1,2) = V.cval(2,2,1,2);      
00781 
00782       (*this)(2,0) = V.cval(1,2,1,1);
00783       (*this)(2,1) = V.cval(1,2,2,2);
00784       (*this)(2,2) = V.cval(1,2,1,2);            
00785 
00786   } else {
00787       if ((numCols != 6) || (numRows != 6)) {      
00788    g3ErrorHandler->warning("Matrix::operator=() - matrix must be %s\n",
00789       "1x1 for tensor of dimension 3\n");      
00790    return *this;
00791       }      
00792       (*this)(0,0) = V.cval(1,1,1,1);
00793       (*this)(0,1) = V.cval(1,1,2,2);
00794       (*this)(0,2) = V.cval(1,1,3,3);      
00795       (*this)(0,3) = V.cval(1,1,1,2);
00796       (*this)(0,4) = V.cval(1,1,1,3);
00797       (*this)(0,5) = V.cval(1,1,2,3);      
00798       
00799       (*this)(1,0) = V.cval(2,2,1,1);
00800       (*this)(1,1) = V.cval(2,2,2,2);
00801       (*this)(1,2) = V.cval(2,2,3,3);      
00802       (*this)(1,3) = V.cval(2,2,1,2);
00803       (*this)(1,4) = V.cval(2,2,1,3);
00804       (*this)(1,5) = V.cval(2,2,2,3);            
00805       
00806       (*this)(2,0) = V.cval(3,3,1,1);
00807       (*this)(2,1) = V.cval(3,3,2,2);
00808       (*this)(2,2) = V.cval(3,3,3,3);      
00809       (*this)(2,3) = V.cval(3,3,1,2);
00810       (*this)(2,4) = V.cval(3,3,1,3);
00811       (*this)(2,5) = V.cval(3,3,2,3);                  
00812       
00813       (*this)(3,0) = V.cval(1,2,1,1);
00814       (*this)(3,1) = V.cval(1,2,2,2);
00815       (*this)(3,2) = V.cval(1,2,3,3);      
00816       (*this)(3,3) = V.cval(1,2,1,2);
00817       (*this)(3,4) = V.cval(1,2,1,3);
00818       (*this)(3,5) = V.cval(1,2,2,3);                        
00819       
00820       (*this)(4,0) = V.cval(1,3,1,1);
00821       (*this)(4,1) = V.cval(1,3,2,2);
00822       (*this)(4,2) = V.cval(1,3,3,3);      
00823       (*this)(4,3) = V.cval(1,3,1,2);
00824       (*this)(4,4) = V.cval(1,3,1,3);
00825       (*this)(4,5) = V.cval(1,3,2,3);                              
00826       
00827       (*this)(5,0) = V.cval(2,3,1,1);
00828       (*this)(5,1) = V.cval(2,3,2,2);
00829       (*this)(5,2) = V.cval(2,3,3,3);      
00830       (*this)(5,3) = V.cval(2,3,1,2);
00831       (*this)(5,4) = V.cval(2,3,1,3);
00832       (*this)(5,5) = V.cval(2,3,2,3);                                    
00833   }
00834   return *this;
00835 }
00836 
00837 
00838 
00839 // virtual Matrix &operator+=(double fact);
00840 // virtual Matrix &operator-=(double fact);
00841 // virtual Matrix &operator*=(double fact);
00842 // virtual Matrix &operator/=(double fact); 
00843 // The above methods all modify the current matrix. If in
00844 // derived matrices data kept in data and of sizeData no redef necessary.
00845 
00846 Matrix &
00847 Matrix::operator+=(double fact)
00848 {
00849   // check if quick return
00850   if (fact == 0.0)
00851     return *this;
00852 
00853   double *dataPtr = data;
00854   for (int i=0; i<dataSize; i++)
00855     *dataPtr++ += fact;
00856   
00857   return *this;
00858 }
00859 
00860 
00861 
00862 
00863 Matrix &
00864 Matrix::operator-=(double fact)
00865 {
00866   // check if quick return
00867   if (fact == 0.0)
00868     return *this;
00869   
00870   double *dataPtr = data;
00871   for (int i=0; i<dataSize; i++)
00872     *dataPtr++ -= fact;
00873 
00874   return *this;
00875 }
00876 
00877 
00878 Matrix &
00879 Matrix::operator*=(double fact)
00880 {
00881   // check if quick return
00882   if (fact == 1.0)
00883     return *this;
00884   
00885   double *dataPtr = data;
00886   for (int i=0; i<dataSize; i++)
00887     *dataPtr++ *= fact;
00888   
00889   return *this;
00890 }
00891 
00892 Matrix &
00893 Matrix::operator/=(double fact)
00894 {
00895     // check if quick return
00896     if (fact == 1.0)
00897  return *this;
00898 
00899     if (fact != 0.0) {
00900       double val = 1.0/fact;
00901 
00902       double *dataPtr = data;
00903       for (int i=0; i<dataSize; i++)
00904  *dataPtr++ *= val;
00905 
00906       return *this;
00907     } else {
00908       // print out the warining message
00909       cerr << "WARNING:Matrix::operator/= - 0 factor specified all values in Matrix set to ";
00910       cerr << MATRIX_VERY_LARGE_VALUE << endl;
00911 
00912       double *dataPtr = data;
00913       for (int i=0; i<dataSize; i++)
00914  *dataPtr++ = MATRIX_VERY_LARGE_VALUE;
00915       
00916       return *this;
00917     }
00918 }
00919 
00920 
00921 //    virtual Matrix operator+(double fact);
00922 //    virtual Matrix operator-(double fact);
00923 //    virtual Matrix operator*(double fact);
00924 //    virtual Matrix operator/(double fact);
00925 // The above methods all return a new full general matrix.
00926 
00927 Matrix
00928 Matrix::operator+(double fact) const
00929 {
00930     Matrix result(*this);
00931     result += fact;
00932     return result;
00933 }
00934 
00935 Matrix
00936 Matrix::operator-(double fact) const
00937 {
00938     Matrix result(*this);
00939     result -= fact;
00940     return result;
00941 }
00942 
00943 Matrix
00944 Matrix::operator*(double fact) const
00945 {
00946     Matrix result(*this);
00947     result *= fact;
00948     return result;
00949 }
00950 
00951 Matrix
00952 Matrix::operator/(double fact) const
00953 {
00954     if (fact == 0.0) {
00955  cerr << "Matrix::operator/(const double &fact): ERROR divide-by-zero\n";
00956  exit(0);
00957     }
00958     Matrix result(*this);
00959     result /= fact;
00960     return result;
00961 }
00962 
00963 
00964 //
00965 // MATRIX_VECTOR OPERATIONS
00966 //
00967 
00968 Vector
00969 Matrix::operator*(const Vector &V) const
00970 {
00971     Vector result(numRows);
00972     
00973     if (V.Size() != numCols) {
00974  cerr << "Matrix::operator*(Vector): incompatable sizes\n";
00975  return result;
00976     } 
00977     
00978     double *dataPtr = data;
00979     for (int i=0; i<numCols; i++)
00980       for (int j=0; j<numRows; j++)
00981  result(j) += *dataPtr++ * V(i);
00982 
00983     /*
00984     cerr << "HELLO: " << V;
00985     for (int i=0; i<numRows; i++) {
00986  double sum = 0.0;
00987  for (int j=0; j<numCols; j++) {
00988      sum += (*this)(i,j) * V(j);
00989      if (i == 9) cerr << "sum: " << sum << " " << (*this)(i,j)*V(j) << " " << V(j) << endl;
00990  }
00991  result(i) += sum;
00992     }
00993     cerr << *this;
00994     cerr << "HELLO result: " << result;    
00995     */
00996 
00997     return result;
00998 }
00999 
01000 Vector
01001 Matrix::operator^(const Vector &V) const
01002 {
01003     Vector result(numCols);
01004     
01005     if (V.Size() != numRows) {
01006       cerr << "Matrix::operator*(Vector): incompatable sizes\n";
01007       return result;
01008     } 
01009 
01010     double *dataPtr = data;
01011     for (int i=0; i<numCols; i++)
01012       for (int j=0; j<numRows; j++)
01013  result(i) += *dataPtr++ * V(j);
01014 
01015     return result;
01016 }
01017 
01018 
01019 //
01020 // MATRIX - MATRIX OPERATIONS
01021 //
01022      
01023 
01024 Matrix
01025 Matrix::operator+(const Matrix &M) const
01026 {
01027     Matrix result(*this);
01028     result.addMatrix(1.0,M,1.0);    
01029     return result;
01030 }
01031      
01032 Matrix
01033 Matrix::operator-(const Matrix &M) const
01034 {
01035     Matrix result(*this);
01036     result.addMatrix(1.0,M,-1.0);    
01037     return result;
01038 }
01039      
01040     
01041 Matrix
01042 Matrix::operator*(const Matrix &M) const
01043 {
01044     Matrix result(numRows,M.numCols);
01045     
01046     if (numCols != M.numRows || result.numRows != numRows) {
01047  cerr << "Matrix::operator*(Matrix): incompatable sizes\n";
01048  return result;
01049     } 
01050 
01051     double *resDataPtr = result.data;     
01052 
01053     int innerDim = numCols;
01054     int nCols = result.numCols;
01055     for (int i=0; i<nCols; i++) {
01056       double *aStartRowDataPtr = data;
01057       double *bStartColDataPtr = &(M.data[i*innerDim]);
01058       for (int j=0; j<numRows; j++) {
01059  double *bDataPtr = bStartColDataPtr;
01060  double *aDataPtr = aStartRowDataPtr +j;     
01061  double sum = 0.0;
01062  for (int k=0; k<innerDim; k++) {
01063    sum += *aDataPtr * *bDataPtr++;
01064    aDataPtr += numRows;
01065  }
01066  *resDataPtr++ = sum;
01067       }
01068     }
01069     return result;
01070 }
01071 
01072 
01073 
01074 // Matrix operator^(const Matrix &M) const
01075 // We overload the * operator to perform matrix^t-matrix multiplication.
01076 // reults = (*this)transposed * M.
01077 
01078 Matrix
01079 Matrix::operator^(const Matrix &M) const
01080 {
01081   Matrix result(numCols,M.numCols);
01082   
01083   if (numRows != M.numRows || result.numRows != numCols) {
01084     cerr << "Matrix::operator*(Matrix): incompatable sizes\n";
01085     return result;
01086   } 
01087 
01088     double *resDataPtr = result.data;     
01089 
01090     int innerDim = numRows;
01091     int nCols = result.numCols;
01092     for (int i=0; i<nCols; i++) {
01093       double *aDataPtr = data;
01094       double *bStartColDataPtr = &(M.data[i*innerDim]);
01095       for (int j=0; j<numCols; j++) {
01096  double *bDataPtr = bStartColDataPtr;
01097  double sum = 0.0;
01098  for (int k=0; k<innerDim; k++) {
01099    sum += *aDataPtr++ * *bDataPtr++;
01100  }
01101  *resDataPtr++ = sum;
01102       }
01103     }
01104 
01105     return result;
01106 }
01107     
01108 
01109 
01110 
01111 Matrix &
01112 Matrix::operator+=(const Matrix &M)
01113 {
01114 #ifdef _G3DEBUG
01115   if (numRows != M.numRows || numCols != M.numCols) {
01116     g3ErrorHandler->warning("Matrix::operator+=(const Matrix &M) - matrices incompatable [%d %d] [%d %d]\n",
01117        numRows, numCols, M.numRows, M.numCols);
01118     return *this;
01119   }
01120 #endif
01121 
01122   double *dataPtr = data;
01123   double *otherData = M.data;
01124   for (int i=0; i<dataSize; i++)
01125     *dataPtr++ += *otherData++;
01126   
01127   return *this;
01128 }
01129 
01130 Matrix &
01131 Matrix::operator-=(const Matrix &M)
01132 {
01133 #ifdef _G3DEBUG
01134   if (numRows != M.numRows || numCols != M.numCols) {
01135     g3ErrorHandler->warning("Matrix::operator-=(const Matrix &M) - matrices incompatable [%d %d] [%d %d]\n",
01136        numRows, numCols, M.numRows, M.numCols);
01137     return *this;
01138   }
01139 #endif
01140 
01141   double *dataPtr = data;
01142   double *otherData = M.data;
01143   for (int i=0; i<dataSize; i++)
01144     *dataPtr++ -= *otherData++;
01145   
01146   return *this;
01147 }
01148 
01149 
01150 //
01151 // Input/Output Methods
01152 //
01153 
01154 void 
01155 Matrix::Output(ostream &s) const
01156 {
01157     for (int i=0; i<noRows(); i++) {
01158  for (int j=0; j<noCols(); j++)
01159      s << (*this)(i,j) << " ";
01160  s << "\n";
01161     }
01162 }
01163 
01164 
01165 void 
01166 Matrix::Input(istream &s)
01167 {
01168     for (int i=0; i<noRows(); i++)
01169  for (int j=0; j<noCols(); j++)
01170      s >> (*this)(i,j);
01171 } 
01172 
01173 
01174 //
01175 // friend stream functions for input and output
01176 //
01177 
01178 ostream &operator<<(ostream &s, const Matrix &V)
01179 {
01180     s << "\n";
01181     V.Output(s);
01182     s << "\n";        
01183     return s;
01184 }
01185 
01186  
01187  
01188 istream &operator>>(istream &s, Matrix &V)
01189 {
01190     V.Input(s);
01191     return s;
01192 }
01193 
01194 
01195 
01196 
01197 
01198 
01199 int
01200 Matrix::Assemble(const Matrix &V, int init_row, int init_col, double fact) 
01201 {
01202   int pos_Rows, pos_Cols;
01203   int res = 0;
01204   
01205   int VnumRows = V.numRows;
01206   int VnumCols = V.numCols;
01207   
01208   int final_row = init_row + VnumRows - 1;
01209   int final_col = init_col + VnumCols - 1;
01210   
01211   if ((init_row >= 0) && (final_row < numRows) && (init_col >= 0) && (final_col < numCols))
01212   {
01213      for (int i=0; i<VnumCols; i++) 
01214      {
01215         pos_Cols = init_col + i;
01216         for (int j=0; j<VnumRows; j++) 
01217         {
01218            pos_Rows = init_row + j;
01219       
01220     (*this)(pos_Rows,pos_Cols) += V(j,i)*fact;
01221         }
01222      }
01223   }  
01224   else 
01225   {
01226      cerr << "WARNING: Matrix::Assemble(const Matrix &V, int init_row, int init_col, double fact): ";
01227      cerr << "position outside bounds \n";
01228      res = -1;
01229   }
01230 
01231   return res;
01232 }
01233 
01234 
01235 
01236 
01237 int
01238 Matrix::AssembleTranspose(const Matrix &V, int init_row, int init_col, double fact) 
01239 {
01240   int pos_Rows, pos_Cols;
01241   int res = 0;
01242   
01243   int VnumRows = V.numRows;
01244   int VnumCols = V.numCols;
01245   
01246   int final_row = init_row + VnumCols - 1;
01247   int final_col = init_col + VnumRows - 1;
01248   
01249   if ((init_row >= 0) && (final_row < numRows) && (init_col >= 0) && (final_col < numCols))
01250   {
01251      for (int i=0; i<VnumRows; i++) 
01252      {
01253         pos_Cols = init_col + i;
01254         for (int j=0; j<VnumCols; j++) 
01255         {
01256            pos_Rows = init_row + j;
01257       
01258     (*this)(pos_Rows,pos_Cols) += V(i,j)*fact;
01259         }
01260      }
01261   }  
01262   else 
01263   {
01264      cerr << "WARNING: Matrix::AssembleTranspose(const Matrix &V, int init_row, int init_col, double fact): ";
01265      cerr << "position outside bounds \n";
01266      res = -1;
01267   }
01268 
01269   return res;
01270 }
01271 
01272 
01273 
01274 
01275 int
01276 Matrix::Extract(const Matrix &V, int init_row, int init_col, double fact) 
01277 {
01278   int pos_Rows, pos_Cols;
01279   int res = 0;
01280   
01281   int VnumRows = V.numRows;
01282   int VnumCols = V.numCols;
01283   
01284   int final_row = init_row + numRows - 1;
01285   int final_col = init_col + numCols - 1;
01286   
01287   if ((init_row >= 0) && (final_row < VnumRows) && (init_col >= 0) && (final_col < VnumCols))
01288   {
01289      for (int i=0; i<numCols; i++) 
01290      {
01291         pos_Cols = init_col + i;
01292         for (int j=0; j<numRows; j++) 
01293         {
01294            pos_Rows = init_row + j;
01295       
01296     (*this)(j,i) = V(pos_Rows,pos_Cols)*fact;
01297         }
01298      }
01299   }  
01300   else 
01301   {
01302      cerr << "WARNING: Matrix::Extract(const Matrix &V, int init_row, int init_col, double fact): ";
01303      cerr << "position outside bounds \n";
01304      res = -1;
01305   }
01306 
01307   return res;
01308 }
01309 
01310 
01311 Matrix operator*(double a, const Matrix &V)
01312 {
01313   return V * a;
01314 }
01315 
01316 
01317 
01318