00001
00002
00003
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
00130 DenseSerialMatrix tmp = *Aptr;
00131
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
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 }