Sandia Home Sandia Home
Main Page | Publications | Downloads | Configuration | Running the Code | Solver Parameters | FAQ | Namespace List | Class Hierarchy | Class List | File List | Namespace Members | Class Members | File Members | Related Pages

APPSPACK_Matrix.cpp

Go to the documentation of this file.
00001 // $Id: APPSPACK_Matrix.cpp,v 1.59 2005/10/16 07:18:36 jgriffi Exp $ 
00002 // $Source: /space/CVS-Acro/acro/packages/appspack/appspack/src/APPSPACK_Matrix.cpp,v $ 
00003 
00004 //@HEADER
00005 // ************************************************************************
00006 // 
00007 //          APPSPACK: Asynchronous Parallel Pattern Search
00008 //                 Copyright (2003) Sandia Corporation
00009 // 
00010 // Under terms of Contract DE-AC04-94AL85000, there is a non-exclusive
00011 // license for use of this work by or on behalf of the U.S. Government.
00012 // 
00013 // This library is free software; you can redistribute it and/or modify
00014 // it under the terms of the GNU Lesser General Public License as
00015 // published by the Free Software Foundation; either version 2.1 of the
00016 // License, or (at your option) any later version.
00017 //  
00018 // This library is distributed in the hope that it will be useful, but
00019 // WITHOUT ANY WARRANTY; without even the implied warranty of
00020 // MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the GNU
00021 // Lesser General Public License for more details.
00022 //                                                                                 
00023 // You should have received a copy of the GNU Lesser General Public
00024 // License along with this library; if not, write to the Free Software
00025 // Foundation, Inc., 59 Temple Place, Suite 330, Boston, MA 02111-1307
00026 // USA.                                                                           .
00027 // 
00028 // Questions? Contact Tammy Kolda (tgkolda@sandia.gov) 
00029 // 
00030 // ************************************************************************
00031 //@HEADER
00032 
00038 #include "APPSPACK_Matrix.hpp"
00039 #include "APPSPACK_Print.hpp"
00040 #include "APPSPACK_Float.hpp"
00041 #include "APPSPACK_LAPACK_Wrappers.hpp"
00042 
00043 APPSPACK::Matrix::Matrix() :
00044   fmatvecSet(false),
00045   fmatvecTSet(false)
00046 {
00047 }
00048 
00049 APPSPACK::Matrix::Matrix(const APPSPACK::Matrix& source, TransposeType ttype) :
00050   fmatvecSet(false),
00051   fmatvecTSet(false)
00052 {
00053   if (ttype == Transpose)
00054     transpose(source);
00055   else
00056     operator=(source);
00057 }
00058 
00059 APPSPACK::Matrix::Matrix(const APPSPACK::Matrix& source, const Vector& s, TransposeType ttype) :
00060   fmatvecSet(false),
00061   fmatvecTSet(false)
00062 {
00063   if (ttype == Transpose)
00064     transpose(source);
00065   else
00066     operator=(source);
00067   scale(s);
00068 }
00069 
00070 APPSPACK::Matrix::Matrix(double** Aptr, int nrows, int ncols) :
00071   fmatvecSet(false),
00072   fmatvecTSet(false)
00073 {
00074   resize(nrows, ncols);
00075   for (int i = 0; i < nrows; i++)
00076     for (int j = 0; j < ncols; j++)
00077       matrix[i][j] = Aptr[i][j];
00078 }
00079 
00080 APPSPACK::Matrix::~Matrix()
00081 {
00082 }
00083 
00084 bool APPSPACK::Matrix::empty() const
00085 {
00086   return (matrix.size() == 0);
00087 }
00088 
00089 int APPSPACK::Matrix::getNrows() const
00090 {
00091   return matrix.size();
00092 }
00093 
00094 int APPSPACK::Matrix::getNcols() const
00095 {
00096   if (matrix.empty())
00097     return 0;
00098 
00099   return matrix[0].size();
00100 }
00101 
00102 const APPSPACK::Vector& APPSPACK::Matrix::getRow(int i) const
00103 {
00104   // Make sure this matrix has at least i rows.
00105   if (getNrows() <= i)
00106   {
00107     cerr << "ERROR: APPSPACK::Matrix::getRow - index out of range." << endl;
00108     throw "APPSPACK Error";
00109   }
00110 
00111   return matrix[i];
00112 }
00113 
00114 void APPSPACK::Matrix::getRowPointers(vector<double *>& Aptr)
00115 { 
00116   for (int i = 0; i < getNrows(); i++)
00117     Aptr.push_back(&matrix[i][0]);
00118   matrixChanged();
00119 }
00120 
00121 void APPSPACK::Matrix::clear()
00122 {
00123   resize(0,0);
00124   matrixChanged();;
00125 }
00126 
00127 void APPSPACK::Matrix::addRow(const APPSPACK::Vector& r)
00128 {
00129   // Make sure length of row r is correct.
00130   if ((!matrix.empty()) && (r.size() != getNcols()))
00131   {
00132     cerr << "ERROR: APPSPACK::Matrix::addRow - incorrect length of input row vector." << endl;
00133     throw "APPSPACK Error";
00134   }
00135   
00136   matrix.push_back(r);
00137   matrixChanged();
00138 }
00139 
00140 void APPSPACK::Matrix::deleteRow(int i)
00141 {
00142   // Make sure this matrix has at least i rows.
00143   if (getNrows() <= i)
00144   {
00145     cerr << "ERROR: APPSPACK::Matrix::getRow - index out of range." << endl;
00146     throw "APPSPACK Error";
00147   }
00148 
00149   matrix.erase(matrix.begin()+i);
00150   matrixChanged();
00151 }
00152 
00153 void APPSPACK::Matrix::addRow(const APPSPACK::Vector& r, double alpha)
00154 {
00155   addRow(r);
00156   matrix[matrix.size()-1].scale(alpha);
00157   matrixChanged();
00158 }
00159 
00160 void APPSPACK::Matrix::addMatrix(const APPSPACK::Matrix& B)
00161 {
00162   for (int i = 0; i < B.getNrows(); i ++)
00163     addRow(B.getRow(i));
00164   matrixChanged();
00165 }
00166 
00167 void APPSPACK::Matrix::addMatrix(const APPSPACK::Matrix& B, double alpha)
00168 {
00169   for (int i = 0; i < B.getNrows(); i ++)
00170     addRow(B.getRow(i), alpha);
00171   matrixChanged();
00172 }
00173 
00174 void APPSPACK::Matrix::addMatrix(const APPSPACK::Matrix& B, const APPSPACK::Vector& s)
00175 {
00176   for (int i = 0; i < B.getNrows(); i ++)
00177   {
00178     addRow(B.getRow(i));
00179     matrix[matrix.size()-1].scale(s);
00180   }
00181   matrixChanged();
00182 }
00183 
00184 void APPSPACK::Matrix::addUniqueRows(const APPSPACK::Matrix& B, double epsilon)
00185 {
00186   // Record number of rows in A now, once we are appending we don't want
00187   // to end up looping over the columns of B also.
00188   int nArows = getNrows();
00189   int nBrows = B.getNrows();
00190   Vector diff(getNcols());
00191   for (int i = 0; i < nBrows; i++)
00192   {
00193     const Vector& bi = B.getRow(i);
00194     bool bi_unique = true;
00195     for (int j = 0; j < nArows; j++)
00196     {
00197       diff = getRow(j);          // diff = aj.
00198       diff -= bi;                // diff = aj - bi.
00199       if (diff.norm() < epsilon) // Is ||aj - bi|| < epsilon?
00200       {
00201         // row is already in A.
00202         bi_unique = false;      
00203         break;
00204       }
00205     }
00206     if (bi_unique)
00207       addRow(bi);
00208   }
00209   
00210   matrixChanged();
00211 }
00212 
00213 void APPSPACK::Matrix::operator=(const APPSPACK::Matrix& m)
00214 {
00215   matrix = m.matrix;
00216   matrixChanged();
00217 }
00218 
00219 void APPSPACK::Matrix::transpose(const Matrix& source)
00220 {
00221   int nrows = source.getNcols();
00222   int ncols = source.getNrows();
00223 
00224   resize(nrows,ncols);
00225 
00226   for (int i = 0; i < nrows; i ++)
00227     for (int j = 0; j < ncols; j ++)
00228       matrix[i][j] = source.matrix[j][i];
00229   matrixChanged();
00230 }
00231 
00232 void APPSPACK::Matrix::scale(const APPSPACK::Matrix& source, const Vector& scaling)
00233 {
00234   // Copy source
00235   matrix = source.matrix; 
00236   
00237   // Scale columns
00238   scale(scaling);
00239   matrixChanged();
00240 }
00241 
00242 void APPSPACK::Matrix::scale(const Vector& scaling)
00243 {
00244   // Scale columns
00245   int nrows = getNrows();
00246   int ncols = getNcols();
00247   for (int i = 0; i < nrows; i++)
00248     for (int j = 0; j < ncols; j++)
00249       matrix[i][j] = matrix[i][j] * scaling[j];
00250   matrixChanged();
00251 }
00252 
00253 void APPSPACK::Matrix::setToIdentity(int n)
00254 {
00255   resize(n,n);
00256   for (int i = 0; i < n; i++)
00257   {
00258     matrix[i].zero();
00259     matrix[i][i] = 1.0;
00260   }
00261   matrixChanged();
00262 }
00263 
00264 void APPSPACK::Matrix::normalize()
00265 {
00266   for (int i = 0; i < getNrows(); i++)
00267   {
00268     double cnorm = matrix[i].norm();
00269     if (cnorm == 0)
00270       deleteRow(i);
00271     else
00272       matrix[i].scale(1/cnorm);
00273   }
00274   matrixChanged();
00275 }
00276 
00277 void  APPSPACK::Matrix::multVec(const APPSPACK::Vector& x, APPSPACK::Vector& y,
00278                                 TransposeType ttype) const
00279 {
00280   // Error check.
00281   if (ttype == Transpose)
00282   {
00283     // Verify length of input vector.
00284     if (x.size() != getNrows())
00285     {
00286       cerr << "ERROR: APPSPACK::Matrix::multVec -  Incorrect length of input vector x." << endl;
00287       throw "APPSPACK Error";
00288     }
00289     
00290     // Verify length of output vector.
00291     if (y.size() != getNcols())
00292     {
00293       cerr << "ERROR: APPSPACK::Matrix::multVec -  Incorrect length of output vector y." << endl;
00294       throw "APPSPACK Error";
00295     }
00296   }
00297   else
00298   {
00299     // Verify length of input vector.
00300     if (x.size() != getNcols())
00301     {
00302       cerr << "ERROR: APPSPACK::Matrix::multVec -  Incorrect length of input vector x." << endl;
00303       throw "APPSPACK Error";
00304     }
00305     
00306     // Verify length of output vector.
00307     if (y.size() != getNrows())
00308     {
00309       cerr << "ERROR: APPSPACK::Matrix::multVec -  Incorrect length of output vector y." << endl;
00310       throw "APPSPACK Error";
00311     }
00312   }
00313 #ifdef HAVE_LAPACK
00314   multVecWithBlas(x,y,ttype);
00315 #else
00316   multVecWithoutBlas(x,y,ttype);
00317 #endif
00318 }
00319 
00320 void APPSPACK::Matrix::svd(APPSPACK::Matrix& U, 
00321                            APPSPACK::Vector &s, 
00322                            APPSPACK::Matrix& VT) const 
00323 {
00324   int m = getNrows();
00325   int n = getNcols();
00326 
00327   // Convert matrix to a FORTRAN vector
00328 
00329   Vector Avec = getMatrixVector();
00330 
00331   // Set up vectors to be filled by the SVD routine
00332 
00333   s.resize(m);                  // Store diag(Sigma)
00334   Vector Uvec(m*m);             // Store U
00335   Vector VTvec(n*n);            // Store V^T
00336 
00337   // Other info
00338 
00339   int info;
00340   char jobu  = 'A';
00341   char jobvt = 'A'; 
00342 
00343   // Note: For call to lapack function dgesvd_, need lwork >=
00344   // max(3*min(m,n)+max(m,n),5*min(m,n)); Works better if slightly
00345   // larger--hence, multiplying quantity by 2.
00346 
00347   int lwork = 2 * max( 3 * min(m,n) + max(m,n), 5 * min(m,n) );
00348   Vector work(lwork);
00349 
00350   DGESVD_F77(&jobu, &jobvt, &m, &n, &Avec[0], &m, &s[0], &Uvec[0], &m,  
00351              &VTvec[0], &n, &work[0], &lwork, &info);
00352 
00353   if (info != 0)
00354   {
00355     cerr << "ERROR: Problem with call to LAPACK function dgesvd.\n";
00356     throw "APPSPACK Error";
00357   }
00358 
00359   U.copyFromFortranVector(Uvec,m,m);
00360   VT.copyFromFortranVector(VTvec,n,n);
00361 }
00362 
00363 void APPSPACK::Matrix::multMat(const APPSPACK::Matrix& B,
00364                                TransposeType ttype)
00365 {
00366   Matrix C;
00367   multMat(B,C,ttype);
00368   matrix = C.matrix;
00369   matrixChanged();
00370 }
00371 
00372 void APPSPACK::Matrix::multMat(const APPSPACK::Matrix& B, APPSPACK::Matrix& C,
00373                                TransposeType ttype) const
00374 {
00375   if (ttype == Transpose)
00376   {
00377     // Peforming mult of form (m x k)*(n x k)' = (m x n).
00378     // Make sure B has correct size.
00379     if (getNcols() != B.getNcols())
00380     {
00381       cerr << "ERROR: Input matrix must have " << getNcols() << " columns." << endl;
00382       throw "APPSPACK Error";
00383     }
00384   }
00385   else
00386   {
00387     // Peforming mult of form (m x k)*(k x n) = (m x n).
00388     // Make sure B has correct size.
00389     if (getNcols() != B.getNrows())
00390     {
00391       cerr << "ERROR: Input matrix must have " << getNrows() << " rows." << endl;
00392       throw "APPSPACK Error";
00393     }
00394   }
00395   
00396 #ifdef HAVE_LAPACK
00397     multMatWithBlas(B,C,ttype);
00398 #else
00399     multMatWithoutBlas(B,C,ttype);
00400 #endif
00401 }
00402 
00403 void APPSPACK::Matrix::nullSpace(APPSPACK::Matrix& ZT, double tol) const
00404 {
00405   // Determine Z such that the columns of Z span the nullspace of this matrix. 
00406 
00407   int m = getNrows();
00408   int n = getNcols();
00409 
00410   if ((m == 0) || (n == 0))
00411   {
00412     cerr << "ERROR: Input matrix must be nonempty." << endl;
00413     throw "APPSPACK Error";
00414   }
00415 
00416   Vector s;
00417   Matrix U;
00418   Matrix VT;
00419 
00420   // Get SVD decomposition
00421   svd(U,s,VT);
00422 
00423   // Determine rank
00424   int rank = s.size();
00425   for (int i = 0; i < s.size(); i++)
00426     if (s[i] < tol)
00427     {
00428       rank = i;
00429       break;
00430     }
00431 
00432   // Form Z' by adding rows of V' that correspond to the zero singular values.
00433   ZT.copySubMatrix(rank,n-rank, VT);
00434 }
00435 
00436 bool APPSPACK::Matrix::getRighInvAndNullBasis(APPSPACK::Matrix& RT, APPSPACK::Matrix& NT,
00437                                             double tol) const
00438 {
00439   RT.clear();
00440   NT.clear();
00441 
00442   int m = getNrows();
00443   int n = getNcols();
00444 
00445   if (m > n)
00446     return false; // Pseudo-inverse does not exist.
00447 
00448   // First compute the singular value decomposition of A.
00449   Matrix U;
00450   Vector s;
00451   Matrix VT;
00452   svd(U,s,VT);
00453   
00454   for (int i = 0; i < s.size(); i++)
00455     if (s[i] < tol)
00456       return false; // Right inverse does not exists.
00457   
00458   // First, let's partition V to be V = [Vr N], i.e., the columns that
00459   // span the range and nullspaces.
00460   Matrix VrT;
00461   VrT.copySubMatrix(0, m, VT);
00462   NT.copySubMatrix(m, n-m, VT);
00463 
00464   // We know that R = Vr*inv(S)*UT => RT = U*inv(S)*VrT.
00465   for (int i = 0; i < s.size(); i++)
00466     s[i] = 1 / s[i];
00467   
00468   RT = U;
00469   RT.scale(s);
00470   RT.multMat(VrT);
00471 
00472   return true;
00473 }
00474 
00475 void APPSPACK::Matrix::pruneDependentRows(APPSPACK::Vector& b, double epsilon)
00476 {
00477   // Perform LQ factorization: A = LQ.
00478   Vector Avec = getMatrixVector(NoTranspose);
00479   int info;
00480   int ncols = getNcols();
00481   int nrows = getNrows();
00482   int lwork = ncols*(ncols+2);
00483   Vector work(lwork, 0.0);
00484   Vector tau(nrows, 0.0);
00485 
00486   DGELQF_F77(&nrows, &ncols, &Avec[0], &nrows, &tau[0], &work[0], &lwork, &info);
00487   
00488   // Delete rows of A and entries of b if necesary.
00489   // Iterate backwards as matrix size may change, i.e. if we were to iterate
00490   // foreward, i for Avec will be different for the i for A and b.  Not to mention nrows.
00491   for (int i = nrows-1; i >= 0; i--)
00492   {
00493     if (fabs(Avec[i+i*nrows]) < epsilon)
00494     {
00495       deleteRow(i);
00496       b.erase(i);
00497     }
00498   }
00499   matrixChanged();
00500 }
00501 
00502 void APPSPACK::Matrix::constrainedLSQR(APPSPACK::Vector& x, const APPSPACK::Vector& b) const
00503 {  
00504   if (empty())
00505     return;
00506 
00507   int info;
00508   int ncols = getNcols();
00509   int nrows = getNrows();
00510   int lwork = ncols*(ncols+2);
00511   Vector work(lwork);
00512   Vector newx(ncols);
00513   Vector bcopy = b;
00514 
00515   Vector Avec = getMatrixVector();
00516 
00517   Vector identity(ncols*ncols, 0.0);
00518   for(int i = 0; i < ncols; i++)
00519     identity[i+i*ncols] = 1;
00520   
00521   DGGLSE_F77(&ncols, &ncols, &nrows, &identity[0], &ncols, &Avec[0], &nrows,
00522              &x[0], &bcopy[0], &newx[0], &work[0], &lwork, &info);
00523 
00524   if (info == 0)
00525     x = newx;
00526 }
00527 
00528 ostream& APPSPACK::Matrix::leftshift(ostream& stream) const
00529 {
00530   stream << endl;
00531   for (int i = 0; i < getNrows(); i++) 
00532   {
00533     for (int j = 0; j < getNcols(); j++) 
00534     {
00535       stream <<  APPSPACK::Print::formatDouble(matrix[i][j]) << " ";
00536     }
00537     stream << endl;
00538   }
00539   return stream;
00540 }
00541 
00542 void APPSPACK::Matrix::matrixChanged()
00543 {
00544   fmatvecTSet = false;
00545   fmatvecSet = false;
00546 }
00547 
00548 const APPSPACK::Vector& APPSPACK::Matrix::getMatrixVector(TransposeType ttype) const
00549 {
00550   if (ttype == Transpose)
00551   {  
00552     if (!fmatvecTSet) // Vectorized matrix not current.  Update.
00553     {
00554       copyToFortranVector(fmatvecT, ttype);
00555       fmatvecTSet = true;
00556     }
00557     return fmatvecT;
00558   }
00559   else
00560   {
00561     if (!fmatvecSet) // Vectorized matrix not current.  Update.
00562     {
00563       copyToFortranVector(fmatvec, ttype);
00564       fmatvecSet = true;
00565     }
00566     return fmatvec;
00567   }
00568 }
00569 
00570 void APPSPACK::Matrix::copyToFortranVector(APPSPACK::Vector& Avec, TransposeType ttype) const
00571 {
00572   // Fortran stores its matrices stored columnwise.
00573   int nrows = getNrows();
00574   int ncols = getNcols();
00575 
00576   Avec.resize(0);
00577   Avec.reserve(nrows * ncols);
00578   
00579   if (ttype == Transpose)
00580   {  
00581     // Store transpose of A. Elements can be added an entire row at a time.
00582     for (int i = 0; i < nrows; i++)
00583       Avec.append(matrix[i]);
00584   }
00585   else
00586   {
00587     // Store A. Elements must be added individually.
00588     for (int j = 0; j < ncols; j++)
00589       for (int i = 0; i < nrows; i ++)
00590         Avec.push_back(matrix[i][j]);
00591   }
00592   
00593 }
00594 
00595 void APPSPACK::Matrix::copyFromFortranVector(const APPSPACK::Vector& Avec, 
00596                                              int nrows, int ncols,
00597                                              TransposeType ttype)
00598 {
00599   // Fortran stores its matrices stored columnwise.
00600 
00601   resize(nrows,ncols);
00602   
00603   if (ttype == Transpose)
00604   {
00605     // Avec contains a transposed version of A. Unpack row-by-row.
00606     int k = 0;
00607     for (int i = 0; i < nrows; i++)
00608       for (int j = 0; j < ncols; j++)
00609         matrix[i][j] = Avec[k++];
00610   }
00611   else
00612   {
00613     // Avec contains A. Unpack column-by-column.
00614     int k = 0;
00615     for (int j = 0; j < ncols; j++)
00616       for (int i = 0; i < nrows; i++)
00617         matrix[i][j] = Avec[k++];
00618   }
00619   
00620   matrixChanged();
00621 }
00622 
00623 // Private
00624 void APPSPACK::Matrix::resize(int nrows, int ncols)
00625 {
00626   matrix.resize(nrows);
00627   for (int i = 0; i < nrows; i++)
00628     matrix[i].resize(ncols);
00629   matrixChanged();
00630 }
00631 
00632 // Private
00633 void APPSPACK::Matrix::copySubMatrix(int istart, int iextent, 
00634                                      const APPSPACK::Matrix& B)
00635 {
00636   if ((istart + iextent) > B.getNrows()) 
00637   {
00638     cerr << "ERROR: Requested submatrix does not exist.";
00639     throw "APPSPACK Error";
00640   }
00641 
00642   clear();
00643   for (int i = 0; i < iextent; i++)
00644     addRow(B.getRow(i+istart));
00645   matrixChanged();
00646 }
00647 
00648 void APPSPACK::Matrix::multMatWithBlas(const APPSPACK::Matrix& B, APPSPACK::Matrix& C,
00649                                        TransposeType ttype) const
00650 {  
00651   int n;  // n denotes number of columns of C (ultimately).
00652           // n is set equal to the number of columns in B if ttype = NoTranspose,
00653           // otherwise, n equals the number of rows in B.
00654  
00655   int ldb; // Leading dimension of B in Fortran lingo.  Essentially,
00656            // ldb equals the number of actual columns in B.  So
00657            // if ttype = NoTranspose, then ldb = n, since B is a
00658            // a (k x n) matrix.  If ttype = Transpose, then ldb = k,
00659            // since in this case B is is a (n x k) matrix.
00660   int m = getNrows();
00661   int k = getNcols();
00662   char transb;
00663   if (ttype == Transpose)
00664   {  
00665     // Peforming mult of form (m x k)*(n x k)' = (m x n).
00666     n      = B.getNrows();
00667     ldb    = k;
00668     transb = 'N';
00669   }
00670   else
00671   {  
00672     // Peforming mult of form (m x k)*(k x n) = (m x n).
00673     n      = B.getNcols();
00674     ldb    = n;
00675     transb = 'T';
00676   }
00677   
00678   // Create memory to hold A*B in vector form.
00679   Vector Cvec(m*n);
00680 
00681 
00682   // We're actually going to compute C = (A')' * (B')' because it's
00683   // more efficient to convert A' and B' to vectors.
00684 
00685   Vector& ATvec = const_cast<Vector&>(getMatrixVector(Transpose));   // Store A' as a vector
00686   Vector& BTvec = const_cast<Vector&>(B.getMatrixVector(Transpose)); // Store B' as a vector
00687   
00688   char transa = 'T';
00689   double one  = 1;
00690   double zero = 0;
00691 
00692   DGEMM_F77(&transa, &transb, &m, &n, &k, &one, &ATvec[0], &k,
00693             &BTvec[0], &ldb , &zero, &Cvec[0], &m);  
00694 
00695   C.copyFromFortranVector(Cvec, m, n);
00696 }
00697 
00698 void APPSPACK::Matrix::multMatWithoutBlas(const APPSPACK::Matrix& B, APPSPACK::Matrix& C,
00699                                           TransposeType ttype) const
00700 {
00701   int m = getNrows();
00702   int k = getNcols();
00703   if (ttype == Transpose)
00704   {
00705     // Peforming mult of form (m x k)*(n x k)' = (m x n).
00706     int n = B.getNrows();
00707     C.resize(m,n);
00708     // C(i,j) = sum_{l=1}^k A(i,l)*B(j,l)
00709     for( int i = 0; i < m; i++)
00710       for( int j = 0; j < n; j++)
00711       {
00712         C.matrix[i][j] = 0;
00713         for( int l = 0; l < k; l++)
00714           C.matrix[i][j] += matrix[i][l] * B.matrix[j][l];
00715       }
00716   }
00717   else
00718   {
00719     // Peforming mult of form (m x k)*(k x n) = (m x n).
00720     int n = B.getNcols();
00721     C.resize(m,n);
00722     // C(i,j) = sum_{l=1}^k A(i,l)*B(l,j)
00723     for( int i = 0; i < m; i++)
00724       for( int j = 0; j < n; j++)
00725       {
00726         C.matrix[i][j] = 0;
00727         for( int l = 0; l < k; l++)
00728           C.matrix[i][j] += matrix[i][l] * B.matrix[l][j];
00729       }
00730   }
00731 }
00732 
00733 void  APPSPACK::Matrix::multVecWithoutBlas(const APPSPACK::Vector& x, APPSPACK::Vector& y,
00734                                            TransposeType ttype) const
00735 {
00736   if (ttype == Transpose)
00737   {
00738     // Compute y = A' * x, i.e. y_j = sum_i A_{ij} x_i
00739     y.zero();
00740     for (int i = 0; i < getNrows(); i++)
00741       for (int j = 0; j < getNcols(); j++)
00742         y[j] += matrix[i][j] * x[i];
00743   }
00744   else
00745   {
00746     // Compute y = A * x, i.e., y_i = sum_j A_{ij} x_j
00747     y.zero();
00748     for (int i = 0; i < getNrows(); i++)
00749       for (int j = 0; j < getNcols(); j++)
00750         y[i] += matrix[i][j] * x[j];
00751   } 
00752 }
00753 
00754 void  APPSPACK::Matrix::multVecWithBlas(const APPSPACK::Vector& x, APPSPACK::Vector& y,
00755                                         TransposeType ttype) const
00756 {
00757   char trans;
00758   if (ttype == Transpose)
00759     trans='T'; // Performing tranpose mult.
00760   else
00761     trans='N'; // Performing normal mult.
00762   int m = getNrows();
00763   int n = getNcols();
00764 
00765   int lda = m;  // Leading dimension.
00766   int incx = 1; // Vectors increment by one.
00767   int incy = 1; // Vectors increment by one.
00768   double beta = 0.0;
00769   double alpha = 1.0;  
00770 
00771   Vector& mcopy = const_cast<Vector&>(getMatrixVector());
00772   Vector& xcopy = const_cast<Vector&>(x);
00773   // Returns y = alpha * mcopy + beta * xcopy.
00774   DGEMV_F77(&trans, &m, &n, &alpha, &mcopy[0], &lda, &xcopy[0], &incx, &beta, &y[0], &incy);  
00775 }
00776 
00777 
00778 ostream& operator<<(ostream& stream, const APPSPACK::Matrix& x)
00779 {
00780   stream << "[";
00781   x.leftshift(stream);
00782   stream << "]";
00783   return stream;
00784 }

 

© Sandia Corporation | Site Contact | Privacy and Security

Generated on Fri Feb 16 10:33:35 2007 for APPSPACK 5.0.1 by doxygen 1.3.9.1 written by Dimitri van Heesch, © 1997-2002