PlayaDenseSerialMatrix.cpp

00001 /* @HEADER@ */
00002 //   
00003  /* @HEADER@ */
00004 
00005 #include "PlayaDenseSerialMatrix.hpp"
00006 #include "PlayaDenseSerialMatrixFactory.hpp"
00007 #include "PlayaSerialVector.hpp"
00008 #include "PlayaVectorSpaceDecl.hpp"  
00009 #include "PlayaVectorDecl.hpp"
00010 #include "PlayaLinearOperatorDecl.hpp"
00011 #include "Teuchos_BLAS.hpp"
00012 
00013 #ifndef HAVE_TEUCHOS_EXPLICIT_INSTANTIATION
00014 #include "PlayaLinearOperatorImpl.hpp"
00015 #include "PlayaVectorImpl.hpp"
00016 #endif
00017 
00018 extern "C"
00019 {
00020 void dgesv_(int *n, int *nrhs, double *a, int* lda, 
00021   int *ipiv, double *b, int *ldb, int *info);
00022 
00023 void dgesvd_( char* jobu, char* jobvt, int* m, int* n, double* a,
00024   int* lda, double* s, double* u, int* ldu, double* vt, int* ldvt,
00025   double* work, int* lwork, int* info );
00026 }
00027 using std::max;
00028 using std::min;
00029 
00030 using namespace Playa;
00031 using namespace Teuchos;
00032 
00033 using std::setw;
00034 
00035 DenseSerialMatrix::DenseSerialMatrix(
00036   const VectorSpace<double>& domain,
00037   const VectorSpace<double>& range)
00038   : LinearOpWithSpaces<double>(domain, range),
00039     nRows_(range.dim()),
00040     nCols_(domain.dim()),
00041     data_(nRows_*nCols_)
00042 {}
00043 
00044 
00045 void DenseSerialMatrix::apply(
00046   Teuchos::ETransp transApplyType,
00047   const Vector<double>& in,
00048   Vector<double> out) const
00049 {
00050   const SerialVector* rvIn = SerialVector::getConcrete(in);
00051   SerialVector* rvOut = SerialVector::getConcrete(out);
00052 
00053   Teuchos::BLAS<int, double> blas;
00054   int lda = numRows();
00055   blas.GEMV(transApplyType, numRows(), numCols(), 1.0, dataPtr(), 
00056     lda, rvIn->dataPtr(), 1, 1.0, rvOut->dataPtr(), 1);
00057 }
00058 
00059 void DenseSerialMatrix::addToRow(int globalRowIndex,
00060   int nElemsToInsert,
00061   const int* globalColumnIndices,
00062   const double* elementValues)
00063 {
00064   int r = globalRowIndex;
00065   for (int k=0; k<nElemsToInsert; k++)
00066   {
00067     int c = globalColumnIndices[k];
00068     double x = elementValues[k];
00069     data_[r + c*numRows()] = x;
00070   }
00071 }
00072 
00073 void DenseSerialMatrix::zero()
00074 {
00075   for (int i=0; i<data_.size(); i++) data_[i] = 0.0;
00076 }
00077 
00078 
00079 void DenseSerialMatrix::print(std::ostream& os) const
00080 {
00081   if (numCols() <= 4)
00082   {
00083     for (int i=0; i<numRows(); i++)
00084     {
00085       for (int j=0; j<numCols(); j++)
00086       {
00087         os << setw(16) << data_[i+numRows()*j];
00088       }
00089       os << std::endl;
00090     }
00091   }
00092   else
00093   {
00094     for (int i=0; i<numRows(); i++)
00095     {
00096       for (int j=0; j<numCols(); j++)
00097       {
00098         os << setw(6) << i << setw(6) << j << setw(20) << data_[i+numRows()*j]
00099            << std::endl;
00100       }
00101     }
00102   }
00103 }
00104 
00105 void DenseSerialMatrix::setRow(int row, const Array<double>& rowVals)
00106 {
00107   TEUCHOS_TEST_FOR_EXCEPT(rowVals.size() != numCols());
00108   TEUCHOS_TEST_FOR_EXCEPT(row < 0);
00109   TEUCHOS_TEST_FOR_EXCEPT(row >= numRows());
00110 
00111   for (int i=0; i<rowVals.size(); i++)
00112   {
00113     data_[row+numRows()*i] = rowVals[i];
00114   }
00115 }
00116 
00117 
00118 namespace Playa
00119 {
00120 
00121 
00122 SolverState<double> denseSolve(const LinearOperator<double>& A,
00123   const Vector<double>& b,
00124   Vector<double>& x)
00125 {
00126   const DenseSerialMatrix* Aptr 
00127     = dynamic_cast<const DenseSerialMatrix*>(A.ptr().get());
00128   TEUCHOS_TEST_FOR_EXCEPT(Aptr==0);
00129   /* make a working copy, because dgesv will overwrite the matrix */
00130   DenseSerialMatrix tmp = *Aptr;
00131   /* Allocate a vector for the solution */
00132   x = b.copy();
00133   SerialVector* xptr 
00134     = dynamic_cast<SerialVector*>(x.ptr().get());
00135   
00136   int N = Aptr->numRows();
00137   int nRHS = 1;
00138   int LDA = N;
00139   Array<int> iPiv(N);
00140   int LDB = N;
00141   int info = 0;
00142   dgesv_(&N, &nRHS, tmp.dataPtr(), &LDA, &(iPiv[0]), xptr->dataPtr(),
00143     &LDB, &info);
00144 
00145   if (info == 0)
00146   {
00147     return SolverState<double>(SolveConverged, "solve OK",
00148       0, 0.0);
00149   }
00150   else 
00151   {
00152     return SolverState<double>(SolveCrashed, "solve crashed with dgesv info="
00153       + Teuchos::toString(info),
00154       0, 0.0);
00155   }
00156 }
00157 
00158 
00159 void denseSVD(const LinearOperator<double>& A,
00160   LinearOperator<double>& U,  
00161   Vector<double>& Sigma,
00162   LinearOperator<double>& Vt)
00163 {
00164   VectorSpace<double> mSpace = A.range();
00165   VectorSpace<double> nSpace = A.domain();
00166 
00167   const DenseSerialMatrix* Aptr 
00168     = dynamic_cast<const DenseSerialMatrix*>(A.ptr().get());
00169   TEUCHOS_TEST_FOR_EXCEPT(Aptr==0);
00170   /* make a working copy, because dgesvd will overwrite the matrix */
00171   DenseSerialMatrix ATmp = *Aptr;
00172 
00173   int M = ATmp.numRows();
00174   int N = ATmp.numCols();
00175   int S = min(M, N);
00176   
00177   VectorSpace<double> sSpace;
00178   if (S==M) sSpace = mSpace;
00179   else sSpace = nSpace;
00180 
00181   Sigma = sSpace.createMember();
00182   SerialVector* sigPtr
00183     = dynamic_cast<SerialVector*>(Sigma.ptr().get());
00184   TEUCHOS_TEST_FOR_EXCEPT(sigPtr==0);
00185 
00186   DenseSerialMatrixFactory umf(sSpace, mSpace);
00187   DenseSerialMatrixFactory vtmf(nSpace, sSpace);
00188   
00189   U = umf.createMatrix();
00190   Vt = vtmf.createMatrix();
00191 
00192   DenseSerialMatrix* UPtr 
00193     = dynamic_cast<DenseSerialMatrix*>(U.ptr().get());
00194   TEUCHOS_TEST_FOR_EXCEPT(UPtr==0);
00195 
00196   DenseSerialMatrix* VtPtr 
00197     = dynamic_cast<DenseSerialMatrix*>(Vt.ptr().get());
00198   TEUCHOS_TEST_FOR_EXCEPT(VtPtr==0);
00199   
00200   double* uData = UPtr->dataPtr();
00201   double* vtData = VtPtr->dataPtr();
00202   double* aData = ATmp.dataPtr();
00203   double* sData = sigPtr->dataPtr();
00204 
00205   char jobu = 'S';
00206   char jobvt = 'S';
00207  
00208   int LDA = M;
00209   int LDU = M;
00210   int LDVT = S;
00211 
00212   int LWORK = max(1, max(3*min(M,N)+max(M,N), 5*min(M,N)));
00213   Array<double> work(LWORK);
00214   
00215   int info = 0;
00216 
00217   dgesvd_(&jobu, &jobvt, &M, &N, aData, &LDA, sData, uData, &LDU, 
00218     vtData, &LDVT, &(work[0]), &LWORK, &info);
00219 
00220   TEUCHOS_TEST_FOR_EXCEPTION(info != 0, std::runtime_error,
00221     "dgesvd failed with error code info=" << info);
00222 
00223   
00224   
00225 }
00226 
00227 }

doxygen