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/07/23 20:41:18 $
00023 // $Source: /usr/local/cvs/OpenSees/SRC/matrix/Vector.cpp,v $
00024                                                                         
00025                                                                         
00026 // File: ~/matrix/Vector.C
00027 //
00028 // Written: fmk 
00029 // Created: 11/96
00030 // Revision: A
00031 //
00032 // Description: This file contains the class implementation for Vector.
00033 //
00034 // What: "@(#) Vector.C, revA"
00035 
00036 #include "Vector.h"
00037 #include "Matrix.h"
00038 #include "ID.h"
00039 
00040 #include <stdlib.h>
00041 #include <math.h>
00042 
00043 double Vector::VECTOR_NOT_VALID_ENTRY =0.0;
00044 
00045 // Vector():
00046 // Standard constructor, sets size = 0;
00047 
00048 Vector::Vector()
00049 : sz(0), theData(0), fromFree(0)
00050 {
00051 
00052 }
00053 
00054 // Vector(int size):
00055 
00056 
00057 
00058 
00059 
00060 
00061 
00062 
00063 
00064 // Constructor used to allocate a vector of size size.
00065 
00066 Vector::Vector(int size)
00067 : sz(size), theData(0), fromFree(0)
00068 {
00069 #ifdef _G3DEBUG
00070   if (sz <= 0) {
00071     g3ErrorHandler->warning("Vector::Vector(int) - size %d specified <= 0\n",size);
00072     sz = 1;
00073   }
00074 #endif
00075 
00076   // get some space for the vector
00077   //  theData = (double *)malloc(size*sizeof(double));
00078   theData = new double [size];
00079 
00080   if (theData == 0) {
00081     g3ErrorHandler->fatal("Vector::Vector(int) - out of memory creating vector of size %d\n",size);
00082     
00083     sz = 0; // set this should fatal error handler not kill process!!
00084   }
00085     
00086   // zero the components
00087   for (int i=0; i<sz; i++)
00088     theData[i] = 0.0;
00089 }
00090 
00091 
00092 // Vector::Vector(double *data, int size)
00093 
00094 Vector::Vector(double *data, int size)
00095 : sz(size),theData(data),fromFree(1)
00096 {
00097 #ifdef _G3DEBUG
00098   if (sz <= 0) {
00099     g3ErrorHandler->warning("Vector::Vector(double *, size) - size %d specified <= 0\n",size);
00100     sz = 0;
00101   }
00102 #endif
00103 }
00104  
00105 
00106 
00107 // Vector(const Vector&):
00108 // Constructor to init a vector from another.
00109 
00110 Vector::Vector(const Vector &other)
00111 : sz(other.sz),theData(0),fromFree(0)
00112 {
00113 #ifdef _G3DEBUG
00114   if (sz < 0) {
00115     g3ErrorHandler->warning("Vector::Vector(int) - size %d specified <= 0\n",sz);
00116     sz = 1;
00117   }
00118 #endif
00119 
00120   //  theData = (double *)malloc(other.sz*sizeof(double));    
00121   theData = new double [other.sz];    
00122   
00123   if (theData == 0) {
00124     g3ErrorHandler->fatal("Vector::Vector(int) - out of memory creating vector of size %d\n",sz);
00125     sz = 0;
00126   }
00127 
00128   // copy the component data
00129   for (int i=0; i<sz; i++)
00130     theData[i] = other.theData[i];
00131 } 
00132 
00133 
00134 
00135 
00136 
00137 
00138 // ~Vector():
00139 //  destructor, deletes the [] data
00140 
00141 Vector::~Vector()
00142 {
00143   if (sz != 0 && fromFree == 0) 
00144     delete [] theData;
00145   //  free((void *)theData);
00146 }
00147 
00148 
00149 int 
00150 Vector::setData(double *newData, int size){
00151   if (sz != 0 && fromFree == 0) 
00152     delete [] theData;      
00153   sz = size;
00154   theData = newData;
00155   fromFree = 1;
00156 
00157   if (sz <= 0) {
00158     g3ErrorHandler->warning("Vector::Vector(double *, size) - size %d specified <= 0\n",size);
00159     sz = 0;
00160   }
00161 
00162   return 0;
00163 }
00164 
00165 
00166 
00167 int 
00168 Vector::resize(int newSize){
00169 
00170   // first check that newSize is valid
00171   if (newSize <= 0) {
00172     g3ErrorHandler->warning("Vector::resize) - size %d specified <= 0\n",newSize); 
00173     return -1;
00174   } 
00175   
00176   // otherwise if newSize is gretaer than oldSize free old space and get new space
00177   else if (newSize > sz) {
00178 
00179     // delete the old array
00180     if (sz != 0 && fromFree == 0) 
00181  delete [] theData;
00182     sz = 0;
00183     fromFree = 0;
00184     
00185     // create new memory
00186     // theData = (double *)malloc(newSize*sizeof(double));    
00187     theData = new double[newSize];
00188     if (theData == 0) {
00189       g3ErrorHandler->fatal("Vector::resize() - out of memory for size %d\n",newSize);
00190       sz = 0;
00191       return -2;
00192     }
00193     sz = newSize;
00194   }  
00195 
00196   // just set the size to be newSize .. penalty of holding onto additional
00197   // memory .. but then save the free() and malloc() calls
00198   else 
00199     sz = newSize;    
00200 
00201   return 0;
00202 }
00203 
00204 
00205 // Assemble(Vector &x, ID &y, double fact ):
00206 // Method to assemble into object the Vector V using the ID l.
00207 // If ID(x) does not exist program writes error message if
00208 // VECTOR_CHECK defined, otherwise ignores it and goes on.
00209 
00210 int 
00211 Vector::Assemble(const Vector &V, const ID &l, double fact )
00212 {
00213   int result = 0;
00214   int pos;
00215   for (int i=0; i<l.Size(); i++) {
00216     pos = l(i);
00217     
00218     if (pos < 0)
00219       ;
00220     else if ((pos < sz) && (i < V.Size()))
00221       // assemble into vector
00222       theData[pos] += V.theData[i] *fact;
00223     else {
00224       result = -1;
00225       if (pos < sz)
00226  g3ErrorHandler->warning("Vector::Assemble() %d out of range [1 %d]\n",
00227     pos, sz-1);
00228       else
00229  g3ErrorHandler->warning("Vector::Assemble() %d out of range [1 %d]\n",
00230     pos, V.Size()-1);
00231     }
00232   }
00233   return result;
00234 }
00235     
00236 
00237 int
00238 Vector::Normalize(void) 
00239 {
00240   double length = 0.0;
00241   for (int i=0; i<sz; i++)
00242     length += theData[i] * theData[i];
00243   length = sqrt(length);
00244   
00245   if (length == 0.0) 
00246     return -1;
00247 
00248   length = 1.0/length;
00249   for (int j=0; j<sz; j++)
00250     theData[j] *= length;  
00251 }
00252 
00253 int
00254 Vector::addVector(double thisFact, const Vector &other, double otherFact )
00255 {
00256   // check if quick return
00257   if (otherFact == 0.0 && thisFact == 1.0)
00258     return 0; 
00259 
00260   // if sizes are compatable add
00261 #ifdef _G3DEBUG
00262   if (sz != other.sz) {
00263     // else sizes are incompatable, do nothing but warning
00264     g3ErrorHandler->warning( "WARNING Vector::addVector() - incompatable Vector sizes\n");
00265     return -1;
00266   }
00267 #endif
00268 
00269 
00270   if (thisFact == 1.0) {
00271 
00272     // want: this += other * otherFact
00273     double *dataPtr = theData;
00274     double *otherDataPtr = other.theData;
00275     if (otherFact == 1.0) { // no point doing a multiplication if otherFact == 1.0
00276       for (int i=0; i<sz; i++) 
00277  *dataPtr++ += *otherDataPtr++;
00278     } else if (otherFact == -1.0) { // no point doing a multiplication if otherFact == 1.0
00279       for (int i=0; i<sz; i++) 
00280  *dataPtr++ -= *otherDataPtr++;
00281     } else 
00282       for (int i=0; i<sz; i++) 
00283  *dataPtr++ += *otherDataPtr++ * otherFact;
00284   } 
00285 
00286   else if (thisFact == 0.0) {
00287 
00288     // want: this = other * otherFact
00289     double *dataPtr = theData;
00290     double *otherDataPtr = other.theData;
00291     if (otherFact == 1.0) { // no point doing a multiplication if otherFact == 1.0
00292       for (int i=0; i<sz; i++) 
00293  *dataPtr++ = *otherDataPtr++;
00294     } else if (otherFact == -1.0) { // no point doing a multiplication if otherFact == 1.0
00295       for (int i=0; i<sz; i++) 
00296  *dataPtr++ = *otherDataPtr++;
00297     } else 
00298       for (int i=0; i<sz; i++) 
00299  *dataPtr++ = *otherDataPtr++ * otherFact;
00300   }
00301 
00302   else {
00303 
00304     // want: this = this * thisFact + other * otherFact
00305     double *dataPtr = theData;
00306     double *otherDataPtr = other.theData;
00307     if (otherFact == 1.0) { // no point doing a multiplication if otherFact == 1.0
00308       for (int i=0; i<sz; i++) {
00309  double value = *dataPtr * thisFact + *otherDataPtr++;
00310  *dataPtr++ = value;
00311       }
00312     } else if (otherFact == -1.0) { // no point doing a multiplication if otherFact == 1.0
00313       for (int i=0; i<sz; i++) {
00314  double value = *dataPtr * thisFact - *otherDataPtr++;
00315  *dataPtr++ = value;
00316       }
00317     } else 
00318       for (int i=0; i<sz; i++) {
00319  double value = *dataPtr * thisFact + *otherDataPtr++ * otherFact;
00320  *dataPtr++ = value;
00321       }
00322   } 
00323 
00324   // successfull
00325   return 0;
00326 }
00327      
00328  
00329 int
00330 Vector::addMatrixVector(double thisFact, const Matrix &m, const Vector &v, double otherFact )
00331 {
00332   // see if quick return
00333   if (thisFact == 1.0 && otherFact == 0.0)
00334     return 0;
00335 
00336   // check the sizes are compatable
00337 #ifdef _G3DEBUG
00338   // check the sizes are compatable
00339   if ((sz != m.noRows()) && (m.noCols() != v.sz)) {
00340     // otherwise incompatable sizes
00341     g3ErrorHandler->warning("Vector::addMatrixVector() %s [%d] = [%d %d] [%d]\n",
00342        "- incompatable sizes",
00343        sz, m.noRows(), m.noCols(), v.sz);
00344     return -1;    
00345   }
00346 #endif
00347 
00348   if (thisFact == 1.0) {
00349 
00350     // want: this += m * v * otherFact
00351     if (otherFact == 1.0) { // no point doing multiplication if otherFact = 1.0
00352       int otherSize = v.sz;
00353       double *matrixDataPtr = m.data;
00354       double *otherDataPtr = v.theData;
00355       for (int i=0; i<otherSize; i++) {
00356  double otherData = *otherDataPtr++;
00357  for (int j=0; j<sz; j++)
00358    theData[j] += *matrixDataPtr++ * otherData;
00359       }
00360     } 
00361     else if (otherFact == -1.0) { // no point doing multiplication if otherFact = -1.0
00362       int otherSize = v.sz;
00363       double *matrixDataPtr = m.data;
00364       double *otherDataPtr = v.theData;
00365       for (int i=0; i<otherSize; i++) {
00366  double otherData = *otherDataPtr++;
00367  for (int j=0; j<sz; j++)
00368    theData[j] -= *matrixDataPtr++ * otherData;
00369       }
00370     } 
00371     else { // have to do the multiplication
00372       int otherSize = v.sz;
00373       double *matrixDataPtr = m.data;
00374       double *otherDataPtr = v.theData;
00375       for (int i=0; i<otherSize; i++) {
00376  double otherData = *otherDataPtr++ * otherFact;
00377  for (int j=0; j<sz; j++)
00378    theData[j] += *matrixDataPtr++ * otherData;
00379       }
00380     }
00381   }
00382 
00383   else if (thisFact == 0.0) {
00384     
00385     // want: this = m * v * otherFact
00386     for (int i=0; i<sz; i++)
00387       theData[i] = 0.0;
00388 
00389     if (otherFact == 1.0) { // no point doing multiplication if otherFact = 1.0
00390       int otherSize = v.sz;
00391       double *matrixDataPtr = m.data;
00392       double *otherDataPtr = v.theData;
00393       for (int i=0; i<otherSize; i++) {
00394  double otherData = *otherDataPtr++;
00395  for (int j=0; j<sz; j++)
00396    theData[j] += *matrixDataPtr++ * otherData;
00397       }
00398     } 
00399     else if (otherFact == -1.0) { // no point doing multiplication if otherFact = -1.0
00400       int otherSize = v.sz;
00401       double *matrixDataPtr = m.data;
00402       double *otherDataPtr = v.theData;
00403       for (int i=0; i<otherSize; i++) {
00404  double otherData = *otherDataPtr++;
00405  for (int j=0; j<sz; j++)
00406    theData[j] -= *matrixDataPtr++ * otherData;
00407       }
00408     } else {
00409       int otherSize = v.sz;
00410       double *matrixDataPtr = m.data;
00411       double *otherDataPtr = v.theData;
00412       for (int i=0; i<otherSize; i++) {
00413  double otherData = *otherDataPtr++ * otherFact;
00414  for (int j=0; j<sz; j++)
00415    theData[j] += *matrixDataPtr++ * otherData;
00416       }
00417     }
00418   }
00419 
00420   else {
00421 
00422     // want: this = this * thisFact + m * v * otherFact
00423     for (int i=0; i<sz; i++)
00424       theData[i] *= thisFact;
00425 
00426     if (otherFact == 1.0) { // no point doing multiplication if otherFact = 1.0
00427       int otherSize = v.sz;
00428       double *matrixDataPtr = m.data;
00429       double *otherDataPtr = v.theData;
00430       for (int i=0; i<otherSize; i++) {
00431  double otherData = *otherDataPtr++;
00432  for (int j=0; j<sz; j++)
00433    theData[j] += *matrixDataPtr++ * otherData;
00434       }
00435     } else if (otherFact == -1.0) { // no point doing multiplication if otherFact = 1.0
00436       int otherSize = v.sz;
00437       double *matrixDataPtr = m.data;
00438       double *otherDataPtr = v.theData;
00439       for (int i=0; i<otherSize; i++) {
00440  double otherData = *otherDataPtr++;
00441  for (int j=0; j<sz; j++)
00442    theData[j] -= *matrixDataPtr++ * otherData;
00443       }
00444     } else {
00445       int otherSize = v.sz;
00446       double *matrixDataPtr = m.data;
00447       double *otherDataPtr = v.theData;
00448       for (int i=0; i<otherSize; i++) {
00449  double otherData = *otherDataPtr++ * otherFact;
00450  for (int j=0; j<sz; j++)
00451    theData[j] += *matrixDataPtr++ * otherData;
00452       }
00453     }
00454   }
00455   
00456   // successfull
00457   return 0;
00458 }
00459 
00460 
00461 
00462 int
00463 Vector::addMatrixTransposeVector(double thisFact, 
00464      const Matrix &m, 
00465      const Vector &v, 
00466      double otherFact )
00467 {
00468   // see if quick return
00469   if (otherFact == 0.0 && thisFact == 1.0)
00470     return 0;
00471 
00472 #ifdef _G3DEBUG
00473   // check the sizes are compatable
00474   if ((sz != m.noRows()) && (m.noRows() != v.sz)) {
00475     // otherwise incompatable sizes
00476     g3ErrorHandler->warning("Vector::addMatrixTransposeVector() %s [%d] = [%d %d] [%d]\n",
00477        "- incompatable sizes",
00478        sz, m.noCols(), m.noRows(), v.sz);
00479     return -1;    
00480   }
00481 #endif
00482 
00483   if (thisFact == 1.0) {
00484 
00485     // want: this += m^t * v * otherFact
00486     if (otherFact == 1.0) { // no point doing multiplication if otherFact = 1.0
00487       int otherSize = v.sz;
00488       double *matrixDataPtr = m.data;
00489       double *otherDataPtrA = v.theData;
00490       for (int i=0; i<sz; i++) {
00491  double *otherDataPtr = otherDataPtrA;
00492  double sum = 0.0;
00493  for (int j=0; j<otherSize; j++)
00494    sum += *matrixDataPtr++ * *otherDataPtr++;
00495  theData[i] += sum;
00496       }
00497     } else if (otherFact == -1.0) { // no point doing multiplication if otherFact = 1.0
00498       int otherSize = v.sz;
00499       double *matrixDataPtr = m.data;
00500       double *otherDataPtrA = v.theData;
00501       for (int i=0; i<sz; i++) {
00502  double *otherDataPtr = otherDataPtrA;
00503  double sum = 0.0;
00504  for (int j=0; j<otherSize; j++)
00505    sum += *matrixDataPtr++ * *otherDataPtr++;
00506  theData[i] -= sum;
00507       }
00508     } else {
00509       int otherSize = v.sz;
00510       double *matrixDataPtr = m.data;
00511       double *otherDataPtrA = v.theData;
00512       for (int i=0; i<sz; i++) {
00513  double *otherDataPtr = otherDataPtrA;
00514  double sum = 0.0;
00515  for (int j=0; j<otherSize; j++)
00516    sum += *matrixDataPtr++ * *otherDataPtr++;
00517  theData[i] += sum * otherFact;
00518       }
00519     }
00520   }
00521 
00522   else if (thisFact == 0.0) {
00523 
00524     // want: this = m^t * v * otherFact
00525     if (otherFact == 1.0) { // no point doing multiplication if otherFact = 1.0
00526       int otherSize = v.sz;
00527       double *matrixDataPtr = m.data;
00528       double *otherDataPtrA = v.theData;
00529       for (int i=0; i<sz; i++) {
00530  double *otherDataPtr = otherDataPtrA;
00531  double sum = 0.0;
00532  for (int j=0; j<otherSize; j++)
00533    sum += *matrixDataPtr++ * *otherDataPtr++;
00534  theData[i] = sum;
00535       }
00536     } else if (otherFact == -1.0) { // no point doing multiplication if otherFact = -1.0
00537       int otherSize = v.sz;
00538       double *matrixDataPtr = m.data;
00539       double *otherDataPtrA = v.theData;
00540       for (int i=0; i<sz; i++) {
00541  double *otherDataPtr = otherDataPtrA;
00542  double sum = 0.0;
00543  for (int j=0; j<otherSize; j++)
00544    sum += *matrixDataPtr++ * *otherDataPtr++;
00545  theData[i] = -sum;
00546       }
00547     } else {
00548       int otherSize = v.sz;
00549       double *matrixDataPtr = m.data;
00550       double *otherDataPtrA = v.theData;
00551       for (int i=0; i<sz; i++) {
00552  double *otherDataPtr = otherDataPtrA;
00553  double sum = 0.0;
00554  for (int j=0; j<otherSize; j++)
00555    sum += *matrixDataPtr++ * *otherDataPtr++;
00556  theData[i] = sum * otherFact;
00557       }
00558     }
00559   } 
00560 
00561   else {
00562 
00563     // want: this = this * thisFact + m^t * v * otherFact
00564     if (otherFact == 1.0) { // no point doing multiplication if otherFact = 1.0
00565       int otherSize = v.sz;
00566       double *matrixDataPtr = m.data;
00567       double *otherDataPtrA = v.theData;
00568       for (int i=0; i<sz; i++) {
00569  double *otherDataPtr = otherDataPtrA;
00570  double sum = 0.0;
00571  for (int j=0; j<otherSize; j++)
00572    sum += *matrixDataPtr++ * *otherDataPtr++;
00573  double value = theData[i] * thisFact + sum;
00574  theData[i] = value;
00575       }
00576     } else if (otherFact == -1.0) { // no point doing multiplication if otherFact = 1.0
00577       int otherSize = v.sz;
00578       double *matrixDataPtr = m.data;
00579       double *otherDataPtrA = v.theData;
00580       for (int i=0; i<sz; i++) {
00581  double *otherDataPtr = otherDataPtrA;
00582  double sum = 0.0;
00583  for (int j=0; j<otherSize; j++)
00584    sum += *matrixDataPtr++ * *otherDataPtr++;
00585  double value = theData[i] * thisFact - sum;
00586  theData[i] = value;
00587       }
00588     } else {
00589       int otherSize = v.sz;
00590       double *matrixDataPtr = m.data;
00591       double *otherDataPtrA = v.theData;
00592       for (int i=0; i<sz; i++) {
00593  double *otherDataPtr = otherDataPtrA;
00594  double sum = 0.0;
00595  for (int j=0; j<otherSize; j++)
00596    sum += *matrixDataPtr++ * *otherDataPtr++;
00597  double value = theData[i] * thisFact + sum * otherFact;
00598  theData[i] = value;
00599       }
00600     }
00601 }
00602 
00603   return 0;
00604 }
00605  
00606  
00607 
00608 
00609 
00610 // double Norm();
00611 // Method to return the norm of  vector. (non-const as may save norm for later)
00612 
00613 double
00614 Vector::Norm(void) const
00615 {
00616   double value = 0;
00617   for (int i=0; i<sz; i++) {
00618     double data = theData[i];
00619     value += data*data;
00620   }
00621   return sqrt(value);
00622 }
00623 
00624 
00625 double &
00626 Vector::operator[](int x) 
00627 {
00628 #ifdef _G3DEBUG
00629   // check if it is inside range [0,sz-1]
00630   if (x < 0 || x >= sz) {
00631       g3ErrorHandler->warning("Vector::() - x %d outside range [0, %d]\n",x,sz-1);
00632       return VECTOR_NOT_VALID_ENTRY;
00633   }
00634 #endif
00635 
00636   return theData[x];
00637 }
00638 
00639 double Vector::operator[](int x) const
00640 {
00641 #ifdef _G3DEBUG
00642   // check if it is inside range [0,sz-1]
00643   if (x < 0 || x >= sz) {
00644     g3ErrorHandler->warning("Vector::() - x %d outside range [0, %d]\n",x,sz-1);
00645     return VECTOR_NOT_VALID_ENTRY;
00646   }
00647 #endif
00648 
00649   return theData[x];
00650 }
00651 
00652 
00653 // operator()(const ID &rows) const
00654 // Method to return a vector whose components are the components of the
00655 // current vector located in positions given by the ID rows.
00656 
00657 
00658 Vector 
00659 Vector::operator()(const ID &rows) const
00660 {
00661   // create a new Vector to be returned
00662   Vector result(rows.Size());
00663 
00664   // check if obtained VEctor of correct size
00665   if (result.Size() != rows.Size()) {
00666     g3ErrorHandler->warning("Vector::()(ID) - new Vector could not be constructed\n");
00667     return result;
00668   }
00669 
00670   // copy the appropraite contents from current to result     
00671   int pos;
00672   for (int i=0; i<rows.Size(); i++) {
00673     pos = rows(i);
00674     if (pos <0 || pos >= sz) {
00675       g3ErrorHandler->warning("Vector::()(ID) - invalid location %d outside range [0,%d]\n",
00676          pos,sz-1);
00677     } else
00678       result(i) = (*this)(pos);
00679   }
00680   return result;
00681 }
00682 
00683 
00684 // Vector &operator=(const Vector  &V):
00685 // the assignment operator, This is assigned to be a copy of V. if sizes
00686 // are not compatable this.theData [] is deleted. The data pointers will not
00687 // point to the same area in mem after the assignment.
00688 //
00689 
00690 Vector &
00691 Vector::operator=(const Vector &V) 
00692 {
00693   // first check we are not trying v = v
00694   if (this != &V) {
00695 
00696    /*
00697 #ifdef _G3DEBUG
00698     // check size compatability, if different warning
00699     if (sz != V.sz) 
00700       g3ErrorHandler->warning("Vector::operator=() - vectors of differing sizes\n");
00701 #endif
00702    */
00703 
00704       if (sz != V.sz)  {
00705 
00706 #ifdef _G3DEBUG
00707    g3ErrorHandler->warning("Vector::operator=() - vectors of differing sizes\n");
00708 #endif
00709 
00710    delete [] this->theData;
00711    this->sz = V.sz;
00712    theData = new double[sz];
00713       }
00714 
00715 
00716       //  copy the data
00717       for (int i=0; i<sz; i++)
00718    theData[i] = V.theData[i];
00719   }
00720 
00721   return *this;
00722 }
00723 
00724 
00725 
00726 
00727 Vector &
00728 Vector::operator=(const Tensor &V) 
00729 {
00730   int rank = V.rank();
00731   if (rank != 2) {
00732       g3ErrorHandler->warning("Vector::operator=() - tensor must be of rank 2\n");
00733       return *this;
00734   }
00735   int dim = V.dim(1);
00736   if (dim != V.dim(2)) {
00737       g3ErrorHandler->warning("Vector::operator=() - tensor must have square dimensions\n");
00738       return *this;
00739   }
00740 
00741   if (dim != 2 || dim != 3 || dim != 1) {
00742       g3ErrorHandler->warning("Vector::operator=() - tensor must be of dimension 2 or 3\n");
00743       return *this;
00744   }      
00745   
00746   if (dim == 1) {
00747       if (sz != 1) {
00748    g3ErrorHandler->warning("Vector::operator=() - Vector size must be 1\n"); 
00749    return *this;
00750       }
00751       theData[0] = V.cval(1,1);
00752   } else if (dim == 2) {
00753       if (sz != 3) {
00754    g3ErrorHandler->warning("Vector::operator=() - Vector size must be 3\n"); 
00755    return *this;
00756       }
00757       theData[0] = V.cval(1,1);
00758       theData[1] = V.cval(2,2);
00759       theData[2] = V.cval(1,2);
00760   } else {
00761       if (sz != 6) {
00762    g3ErrorHandler->warning("Vector::operator=() - Vector size must be 6\n"); 
00763    return *this;
00764       }      
00765       theData[0] = V.cval(1,1);
00766       theData[1] = V.cval(2,2);
00767       theData[2] = V.cval(3,3);
00768       theData[3] = V.cval(1,2);
00769       theData[4] = V.cval(1,3);
00770       theData[5] = V.cval(2,3);
00771   }
00772   return *this;
00773 }
00774 
00775 
00776 
00777 
00778 // Vector &operator+=(double fact):
00779 // The += operator adds fact to each element of the vector, data[i] = data[i]+fact.
00780 
00781 Vector &Vector::operator+=(double fact)
00782 {
00783   if (fact != 0.0)
00784     for (int i=0; i<sz; i++)
00785       theData[i] += fact;
00786   return *this;    
00787 }
00788 
00789 
00790 
00791 // Vector &operator-=(double fact)
00792 // The -= operator subtracts fact from each element of the vector, data[i] = data[i]-fact.
00793 
00794 Vector &Vector::operator-=(double fact)
00795 {
00796   if (fact != 0.0)
00797     for (int i=0; i<sz; i++)
00798       theData[i] -= fact;
00799   return *this;    
00800 }
00801 
00802 
00803 
00804 // Vector &operator*=(double fact):
00805 // The -= operator subtracts fact from each element of the vector, data[i] = data[i]-fact.
00806 
00807 Vector &Vector::operator*=(double fact)
00808 {
00809   for (int i=0; i<sz; i++)
00810     theData[i] *= fact;
00811   return *this;
00812 }
00813 
00814 
00815 
00816 // Vector &operator/=(double fact):
00817 // The /= operator divides each element of the vector by fact, theData[i] = theData[i]/fact.
00818 // Program exits if divide-by-zero would occur with warning message.
00819 
00820 Vector &Vector::operator/=(double fact)
00821 {
00822   if (fact == 0.0) { // instead of divide-by-zero error set to VECTOR_VERY_LARGE_VALUE
00823     for (int i=0; i<sz; i++)
00824       theData[i] = VECTOR_VERY_LARGE_VALUE;
00825   } else {
00826     for (int i=0; i<sz; i++)
00827       theData[i] /= fact;
00828   }
00829   
00830   return *this;
00831 }
00832 
00833 
00834 
00835 
00836 // Vector operator+(double fact):
00837 // The + operator returns a Vector of the same size as current, whose components
00838 // are return(i) = theData[i]+fact;
00839 
00840 Vector 
00841 Vector::operator+(double fact) const
00842 {
00843   Vector result(*this);
00844   if (result.Size() != sz) 
00845     g3ErrorHandler->warning("Vector::operator+(double) - ran out of memory for new Vector\n");
00846 
00847   result += fact;
00848   return result;
00849 }
00850 
00851 
00852 
00853 // Vector operator-(double fact):
00854 // The + operator returns a Vector of the same size as current, whose components
00855 // are return(i) = theData[i]-fact;
00856 
00857 Vector 
00858 Vector::operator-(double fact) const
00859 {
00860     Vector result(*this);
00861     if (result.Size() != sz) 
00862       g3ErrorHandler->warning("Vector::operator-(double) - ran out of memory for new Vector\n");
00863 
00864     result -= fact;
00865     return result;
00866 }
00867 
00868 
00869 
00870 // Vector operator*(double fact):
00871 // The + operator returns a Vector of the same size as current, whose components
00872 // are return(i) = theData[i]*fact;
00873 
00874 Vector 
00875 Vector::operator*(double fact) const
00876 {
00877     Vector result(*this);
00878     if (result.Size() != sz) 
00879       g3ErrorHandler->warning("Vector::operator-(double) - ran out of memory for new Vector\n");
00880 
00881     result *= fact;
00882     return result;
00883 }
00884 
00885 
00886 // Vector operator/(double fact):
00887 // The + operator returns a Vector of the same size as current, whose components
00888 // are return(i) = theData[i]/fact; Exits if divide-by-zero error.
00889 
00890 Vector 
00891 Vector::operator/(double fact) const
00892 {
00893     if (fact == 0.0) 
00894       g3ErrorHandler->warning("Vector::operator/(double fact) - divide-by-zero error coming\n");
00895 
00896     Vector result(*this);
00897     if (result.Size() != sz) 
00898       g3ErrorHandler->warning("Vector::operator-(double) - ran out of memory for new Vector\n");
00899 
00900     result /= fact;
00901     return result;
00902 }
00903 
00904 
00905 
00906 // Vector &operator+=(const Vector &V):
00907 // The += operator adds V's data to data, data[i]+=V(i). A check to see if
00908 // vectors are of same size is performed if VECTOR_CHECK is defined.
00909 
00910 Vector &
00911 Vector::operator+=(const Vector &other)
00912 {
00913 #ifdef _G3DEBUG
00914   if (sz != other.sz) {
00915     g3ErrorHandler->warning("WARNING Vector::operator+=(Vector):Vectors not of same sizes: %d %d\n",
00916         sz, other.sz);
00917     return *this;
00918   }    
00919 #endif
00920 
00921   for (int i=0; i<sz; i++)
00922     theData[i] += other.theData[i];
00923   return *this;     
00924 }
00925 
00926 
00927 
00928 // Vector &operator-=(const Vector &V):
00929 // The -= operator subtracts V's data from  data, data[i]+=V(i). A check 
00930 //    to see if vectors are of same size is performed if VECTOR_CHECK is defined.
00931 
00932 Vector &
00933 Vector::operator-=(const Vector &other)
00934 {
00935 #ifdef _G3DEBUG
00936   if (sz != other.sz) {
00937     g3ErrorHandler->warning("WARNING Vector::operator-=(Vector): Vectors not of same sizes: %d %d\n",
00938        sz, other.sz);
00939     return *this;
00940   }
00941 #endif
00942   
00943   for (int i=0; i<sz; i++)
00944     theData[i] -= other.theData[i];
00945   return *this;    
00946 }
00947 
00948 
00949 
00950 // Vector operator+(const Vector &V):
00951 // The + operator checks the two vectors are of the same size if VECTOR_CHECK is defined.
00952 //  Then returns a Vector whose components are the vector sum of current and V's data.
00953 
00954 Vector 
00955 Vector::operator+(const Vector &b) const
00956 {
00957 #ifdef _G3DEBUG
00958   if (sz != b.sz) {
00959     g3ErrorHandler->warning( "Vector::operator+(Vector): Vectors not of same size, sizes: %d %d",
00960         sz, b.sz);
00961         
00962     return *this;
00963   }
00964 #endif
00965 
00966     Vector result(*this);
00967 
00968     // check new Vector of correct size
00969   if (result.Size() != sz) {
00970     g3ErrorHandler->warning("Vector::operator-(Vector): new Vector not of correct size ");
00971     return result;
00972   }
00973   result += b;
00974   return result;
00975 }
00976 
00977 
00978 // Vector operator-(const Vector &V):
00979 // The - operator checks the two vectors are of the same size and then returns a Vector
00980 // whose components are the vector difference of current and V's data.
00981 
00982 Vector 
00983 Vector::operator-(const Vector &b) const
00984 {
00985 #ifdef _G3DEBUG
00986   if (sz != b.sz) {
00987     g3ErrorHandler->warning("WARNING Vector::operator-(Vector): Vectors not of same sizes: %d %d\n",
00988        sz, b.sz);
00989     return *this;
00990   }
00991 #endif
00992 
00993   Vector result(*this);
00994 
00995   // check new Vector of correct size
00996   if (result.Size() != sz) {
00997     g3ErrorHandler->warning("Vector::operator-(Vector): new Vector not of correct size ");
00998     return result;
00999   }
01000 
01001   result -= b;
01002   return result;
01003 }
01004 
01005 
01006 
01007 // double operator^(const Vector &V) const;
01008 // Method to perform (Vector)transposed * vector.
01009 double
01010 Vector::operator^(const Vector &V) const
01011 {
01012 #ifdef _G3DEBUG
01013   if (sz != V.sz) {
01014     g3ErrorHandler->warning( "WARNING Vector::operator-(Vector): Vectors not of same sizes: %d %d\n",
01015        sz, V.sz);
01016     return 0.0;
01017   }
01018 #endif
01019 
01020   double result = 0.0;
01021   double *dataThis = theData;
01022   double *dataV = V.theData;
01023     for (int i=0; i<sz; i++)
01024       result += *dataThis++ * *dataV++;
01025 
01026     return result;
01027 }
01028 
01029 
01030 // Vector operator/(const Matrix &M) const;    
01031 // Method to return inv(M)*this
01032 
01033 Vector
01034 Vector::operator/(const Matrix &M) const
01035 {
01036   Vector res(M.noRows());
01037     
01038   if (M.noRows() != M.noCols()) { // if not square do least squares solution
01039     Matrix A(M^M);
01040     A.Solve(*this, res);    
01041   }
01042   else {
01043     M.Solve(*this, res);
01044   }
01045   return res;
01046 }
01047     
01048  
01049 
01050  
01051 
01052 
01053 
01054 // friend ostream &operator<<(ostream &s, const Vector &V)
01055 // A function is defined to allow user to print the vectors using ostreams.
01056 
01057 ostream &operator<<(ostream &s, const Vector &V)
01058 {
01059   for (int i=0; i<V.Size(); i++) 
01060       s << V(i) << " ";
01061 
01062   return s << "\n";
01063 }
01064 
01065 // friend istream &operator>>(istream &s, Vector &V)
01066 // A function is defined to allow user to input the data into a Vector which has already
01067 // been constructed with data, i.e. Vector(int) or Vector(const Vector &) constructors.
01068 
01069 istream &operator>>(istream &s, Vector &V)
01070 {
01071   for (int i=0; i<V.Size(); i++) 
01072     s >> V(i);
01073   
01074     return s;
01075 }
01076 
01077 
01078 
01079 Vector operator*(double a, const Vector &V)
01080 {
01081   return V * a;
01082 }
01083 
01084 
01085 int
01086 Vector::Assemble(const Vector &V, int init_pos, double fact) 
01087 {
01088   int res = 0;
01089   int cur_pos   = init_pos;  
01090   int final_pos = init_pos + V.sz - 1;
01091   
01092   if ((init_pos >= 0) && (final_pos < sz))
01093   {
01094      for (int j=0; j<V.sz; j++) 
01095         (*this)(cur_pos++) += V(j)*fact;
01096   }
01097   else 
01098   {
01099      cerr << "WARNING: Vector::Assemble(const Vector &V, int init_pos, double fact): ";
01100      cerr << "position outside bounds \n";
01101      res = -1;
01102   }
01103 
01104   return res;
01105 }
01106 
01107 
01108 
01109 
01110 int
01111 Vector::Extract(const Vector &V, int init_pos, double fact) 
01112 {
01113   int res = 0;
01114   int cur_pos   = init_pos;  
01115   int final_pos = init_pos + sz - 1;
01116   
01117   if ((init_pos >= 0) && (final_pos < V.sz))
01118   {
01119      for (int j=0; j<sz; j++) 
01120         (*this)(j) = V(cur_pos++)*fact;
01121   }
01122   else 
01123   {
01124      cerr << "WARNING: Vector::Assemble(const Vector &V, int init_pos, double fact): ";
01125      cerr << "position outside bounds \n";
01126      res = -1;
01127   }
01128 
01129   return res;
01130 }