DenseLinAlgLAPack.cpp

Go to the documentation of this file.
00001 // @HEADER
00002 // ***********************************************************************
00003 // 
00004 // Moocho: Multi-functional Object-Oriented arCHitecture for Optimization
00005 //                  Copyright (2003) Sandia Corporation
00006 // 
00007 // Under terms of Contract DE-AC04-94AL85000, there is a non-exclusive
00008 // license for use of this work by or on behalf of the U.S. Government.
00009 // 
00010 // This library is free software; you can redistribute it and/or modify
00011 // it under the terms of the GNU Lesser General Public License as
00012 // published by the Free Software Foundation; either version 2.1 of the
00013 // License, or (at your option) any later version.
00014 //  
00015 // This library is distributed in the hope that it will be useful, but
00016 // WITHOUT ANY WARRANTY; without even the implied warranty of
00017 // MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the GNU
00018 // Lesser General Public License for more details.
00019 //  
00020 // You should have received a copy of the GNU Lesser General Public
00021 // License along with this library; if not, write to the Free Software
00022 // Foundation, Inc., 59 Temple Place, Suite 330, Boston, MA 02111-1307
00023 // USA
00024 // Questions? Contact Roscoe A. Bartlett (rabartl@sandia.gov) 
00025 // 
00026 // ***********************************************************************
00027 // @HEADER
00028 
00029 #include "DenseLinAlgLAPack.hpp"
00030 #include "DenseLinAlgPack_LAPACK_Cpp.hpp"
00031 #include "DenseLinAlgPack_DMatrixAsTriSym.hpp"
00032 #include "Teuchos_TestForException.hpp"
00033 
00034 namespace {
00035 template< class T >
00036 inline
00037 T my_min( const T& v1, const T& v2 ) { return v1 < v2 ? v1 : v2; }
00038 } // end namespace
00039 
00040 void DenseLinAlgLAPack::potrf( DMatrixSliceTriEle* A )
00041 {
00042   FortranTypes::f_int info;
00043   LAPACK_Cpp::potrf( A->uplo(), A->rows(), A->gms().col_ptr(1), A->gms().max_rows(), &info );
00044   if( info != 0 ) {
00045     TEST_FOR_EXCEPTION(
00046       info < 0, std::invalid_argument
00047       ,"potrf(...): Error, Invalid argument "
00048       << -info << " sent to LAPACK function xPOTRF(...)" );
00049     // info > 0
00050     TEST_FOR_EXCEPTION(
00051       true, FactorizationException
00052       ,"potrf(...): Error, Minor of order "
00053       << info << " is not positive definite, the factorization "
00054       "could not be completed" );
00055   }
00056   // If you get here all went okay and A has been factorized and now
00057   // hold the cholesky factor of input A.
00058 }
00059 
00060 void DenseLinAlgLAPack::geqrf( DMatrixSlice* A, DVectorSlice* tau, DVectorSlice* work )
00061 {
00062   FortranTypes::f_int info;
00063   if( tau->dim() != my_min( A->rows(), A->cols() ) ) {
00064     TEST_FOR_EXCEPTION(
00065       true, std::invalid_argument, "geqrf(...): Error, tau is not sized correctly!" );
00066   }
00067   LAPACK_Cpp::geqrf( A->rows(), A->cols(),A->col_ptr(1), A->max_rows()
00068     , tau->raw_ptr(), work->raw_ptr(), work->dim(), &info  );
00069   if( info != 0 ) {
00070     TEST_FOR_EXCEPTION(
00071       info < 0, std::invalid_argument
00072       ,"geqrf(...): Error, Invalid argument "
00073       << -info << " sent to LAPACK function xGEQRF(...)" );
00074   }
00075   // If you get here A and tau contain the QR factors of input A.
00076 }
00077 
00078 void DenseLinAlgLAPack::ormrq(
00079   BLAS_Cpp::Side side, BLAS_Cpp::Transp trans
00080   ,const DMatrixSlice& A, const DVectorSlice& tau
00081   ,DMatrixSlice* C, DVectorSlice* work
00082   )
00083 {
00084   FortranTypes::f_int info;
00085   LAPACK_Cpp::ormqr( side, trans, C->rows(), C->cols()
00086     , tau.dim(), A.col_ptr(1), A.max_rows()
00087     , tau.raw_ptr(), C->col_ptr(1), C->max_rows()
00088     , work->raw_ptr(), work->dim(), &info );
00089   if( info != 0 ) {
00090     TEST_FOR_EXCEPTION(
00091       info < 0, std::invalid_argument
00092       ,"ormrq(...): Error, Invalid argument "
00093       << -info << " sent to LAPACK function xORMRQ(...)" );
00094   }
00095   // If you get here C contains the desired matrix-matrix multiplication.
00096 }
00097 
00098 void DenseLinAlgLAPack::sytrf(
00099   DMatrixSliceTriEle* A, FortranTypes::f_int ipiv[]
00100   ,DVectorSlice* work
00101   )
00102 {
00103   FortranTypes::f_int info;
00104   LAPACK_Cpp::sytrf( A->uplo(), A->rows(), A->gms().col_ptr(1)
00105     , A->gms().max_rows(), ipiv, work->raw_ptr(), work->dim()
00106     , &info );
00107   TEST_FOR_EXCEPTION(
00108     info < 0, std::invalid_argument
00109     ,"sytrf(...): Error, Invalid argument "
00110     << -info << " sent to LAPACK function xSYTRF(...)" );
00111   TEST_FOR_EXCEPTION(
00112     info > 0, FactorizationException
00113     ,"sytrf(...): Error, xSYTRF(...) indicates a singular matrix, "
00114     << "D("<<info<<","<<info<<") is zero." );
00115   // If we get here A and ipiv contain the factorization.
00116 }
00117 
00118 void DenseLinAlgLAPack::sytrs(
00119   const DMatrixSliceTriEle& A, FortranTypes::f_int ipiv[]
00120   ,DMatrixSlice* B, DVectorSlice* work
00121   )
00122 {
00123   TEST_FOR_EXCEPTION(
00124     (A.rows() != B->rows()), std::invalid_argument
00125     ,"sytrs(...) : Error, The number of rows in A and B must match."
00126     );
00127   FortranTypes::f_int info;
00128   LAPACK_Cpp::sytrs(  A.uplo(), A.rows(), B->cols(), A.gms().col_ptr(1)
00129     , A.gms().max_rows(), ipiv, B->col_ptr(1), B->max_rows()
00130     , &info );
00131   TEST_FOR_EXCEPTION(
00132     info < 0, std::invalid_argument
00133     ,"sytrs(...): Error, Invalid argument "
00134     << -info << " sent to LAPACK function xSYTRS(...)"
00135     );
00136 }
00137 
00138 void DenseLinAlgLAPack::getrf(
00139   DMatrixSlice* A, FortranTypes::f_int ipiv[], FortranTypes::f_int* rank
00140   )
00141 {
00142   FortranTypes::f_int info;
00143   LAPACK_Cpp::getrf(
00144     A->rows(), A->cols(), A->col_ptr(1), A->max_rows()
00145     ,ipiv, &info
00146     );
00147   *rank = my_min( A->rows(), A->cols() );
00148   TEST_FOR_EXCEPTION(
00149     info < 0, std::invalid_argument
00150     ,"getrf(...): Error, Invalid argument "
00151     << -info << " sent to LAPACK function xGETRF(...)" );
00152   if(info > 0)
00153     *rank = info - 1;
00154 }
00155 
00156 void DenseLinAlgLAPack::getrs(
00157   const DMatrixSlice& LU, const FortranTypes::f_int ipiv[], BLAS_Cpp::Transp transp
00158   ,DMatrixSlice* B
00159   )
00160 {
00161   TEST_FOR_EXCEPTION(
00162     (LU.rows() != LU.cols() || LU.rows() != B->rows() ), std::invalid_argument
00163     ,"getrs(...) : Error, A must be square and the number of rows in A and B must match."
00164     );
00165   FortranTypes::f_int info;
00166   LAPACK_Cpp::getrs(
00167     transp, LU.rows(), B->cols(), LU.col_ptr(1), LU.max_rows(), ipiv
00168     ,B->col_ptr(1), B->max_rows(), &info
00169     );
00170   TEST_FOR_EXCEPTION(
00171     info < 0, std::invalid_argument
00172     ,"getrs(...): Error, Invalid argument "
00173     << -info << " sent to LAPACK function xGETRS(...)" );
00174   // If we get here B is the solution to the linear system.
00175 }

Generated on Wed May 12 21:52:29 2010 for MOOCHO (Single Doxygen Collection) by  doxygen 1.4.7