MOOCHO (Single Doxygen Collection) Version of the Day
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 // Redistribution and use in source and binary forms, with or without
00011 // modification, are permitted provided that the following conditions are
00012 // met:
00013 //
00014 // 1. Redistributions of source code must retain the above copyright
00015 // notice, this list of conditions and the following disclaimer.
00016 //
00017 // 2. Redistributions in binary form must reproduce the above copyright
00018 // notice, this list of conditions and the following disclaimer in the
00019 // documentation and/or other materials provided with the distribution.
00020 //
00021 // 3. Neither the name of the Corporation nor the names of the
00022 // contributors may be used to endorse or promote products derived from
00023 // this software without specific prior written permission.
00024 //
00025 // THIS SOFTWARE IS PROVIDED BY SANDIA CORPORATION "AS IS" AND ANY
00026 // EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED TO, THE
00027 // IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR
00028 // PURPOSE ARE DISCLAIMED. IN NO EVENT SHALL SANDIA CORPORATION OR THE
00029 // CONTRIBUTORS BE LIABLE FOR ANY DIRECT, INDIRECT, INCIDENTAL, SPECIAL,
00030 // EXEMPLARY, OR CONSEQUENTIAL DAMAGES (INCLUDING, BUT NOT LIMITED TO,
00031 // PROCUREMENT OF SUBSTITUTE GOODS OR SERVICES; LOSS OF USE, DATA, OR
00032 // PROFITS; OR BUSINESS INTERRUPTION) HOWEVER CAUSED AND ON ANY THEORY OF
00033 // LIABILITY, WHETHER IN CONTRACT, STRICT LIABILITY, OR TORT (INCLUDING
00034 // NEGLIGENCE OR OTHERWISE) ARISING IN ANY WAY OUT OF THE USE OF THIS
00035 // SOFTWARE, EVEN IF ADVISED OF THE POSSIBILITY OF SUCH DAMAGE.
00036 //
00037 // Questions? Contact Roscoe A. Bartlett (rabartl@sandia.gov) 
00038 // 
00039 // ***********************************************************************
00040 // @HEADER
00041 
00042 #include "DenseLinAlgLAPack.hpp"
00043 #include "DenseLinAlgPack_LAPACK_Cpp.hpp"
00044 #include "DenseLinAlgPack_DMatrixAsTriSym.hpp"
00045 #include "Teuchos_Assert.hpp"
00046 
00047 namespace {
00048 template< class T >
00049 inline
00050 T my_min( const T& v1, const T& v2 ) { return v1 < v2 ? v1 : v2; }
00051 } // end namespace
00052 
00053 void DenseLinAlgLAPack::potrf( DMatrixSliceTriEle* A )
00054 {
00055   FortranTypes::f_int info;
00056   LAPACK_Cpp::potrf( A->uplo(), A->rows(), A->gms().col_ptr(1), A->gms().max_rows(), &info );
00057   if( info != 0 ) {
00058     TEUCHOS_TEST_FOR_EXCEPTION(
00059       info < 0, std::invalid_argument
00060       ,"potrf(...): Error, Invalid argument "
00061       << -info << " sent to LAPACK function xPOTRF(...)" );
00062     // info > 0
00063     TEUCHOS_TEST_FOR_EXCEPTION(
00064       true, FactorizationException
00065       ,"potrf(...): Error, Minor of order "
00066       << info << " is not positive definite, the factorization "
00067       "could not be completed" );
00068   }
00069   // If you get here all went okay and A has been factorized and now
00070   // hold the cholesky factor of input A.
00071 }
00072 
00073 void DenseLinAlgLAPack::geqrf( DMatrixSlice* A, DVectorSlice* tau, DVectorSlice* work )
00074 {
00075   FortranTypes::f_int info;
00076   if( tau->dim() != my_min( A->rows(), A->cols() ) ) {
00077     TEUCHOS_TEST_FOR_EXCEPTION(
00078       true, std::invalid_argument, "geqrf(...): Error, tau is not sized correctly!" );
00079   }
00080   LAPACK_Cpp::geqrf( A->rows(), A->cols(),A->col_ptr(1), A->max_rows()
00081     , tau->raw_ptr(), work->raw_ptr(), work->dim(), &info  );
00082   if( info != 0 ) {
00083     TEUCHOS_TEST_FOR_EXCEPTION(
00084       info < 0, std::invalid_argument
00085       ,"geqrf(...): Error, Invalid argument "
00086       << -info << " sent to LAPACK function xGEQRF(...)" );
00087   }
00088   // If you get here A and tau contain the QR factors of input A.
00089 }
00090 
00091 void DenseLinAlgLAPack::ormrq(
00092   BLAS_Cpp::Side side, BLAS_Cpp::Transp trans
00093   ,const DMatrixSlice& A, const DVectorSlice& tau
00094   ,DMatrixSlice* C, DVectorSlice* work
00095   )
00096 {
00097   FortranTypes::f_int info;
00098   LAPACK_Cpp::ormqr( side, trans, C->rows(), C->cols()
00099     , tau.dim(), A.col_ptr(1), A.max_rows()
00100     , tau.raw_ptr(), C->col_ptr(1), C->max_rows()
00101     , work->raw_ptr(), work->dim(), &info );
00102   if( info != 0 ) {
00103     TEUCHOS_TEST_FOR_EXCEPTION(
00104       info < 0, std::invalid_argument
00105       ,"ormrq(...): Error, Invalid argument "
00106       << -info << " sent to LAPACK function xORMRQ(...)" );
00107   }
00108   // If you get here C contains the desired matrix-matrix multiplication.
00109 }
00110 
00111 void DenseLinAlgLAPack::sytrf(
00112   DMatrixSliceTriEle* A, FortranTypes::f_int ipiv[]
00113   ,DVectorSlice* work
00114   )
00115 {
00116   FortranTypes::f_int info;
00117   LAPACK_Cpp::sytrf( A->uplo(), A->rows(), A->gms().col_ptr(1)
00118     , A->gms().max_rows(), ipiv, work->raw_ptr(), work->dim()
00119     , &info );
00120   TEUCHOS_TEST_FOR_EXCEPTION(
00121     info < 0, std::invalid_argument
00122     ,"sytrf(...): Error, Invalid argument "
00123     << -info << " sent to LAPACK function xSYTRF(...)" );
00124   TEUCHOS_TEST_FOR_EXCEPTION(
00125     info > 0, FactorizationException
00126     ,"sytrf(...): Error, xSYTRF(...) indicates a singular matrix, "
00127     << "D("<<info<<","<<info<<") is zero." );
00128   // If we get here A and ipiv contain the factorization.
00129 }
00130 
00131 void DenseLinAlgLAPack::sytrs(
00132   const DMatrixSliceTriEle& A, FortranTypes::f_int ipiv[]
00133   ,DMatrixSlice* B, DVectorSlice* work
00134   )
00135 {
00136   TEUCHOS_TEST_FOR_EXCEPTION(
00137     (A.rows() != B->rows()), std::invalid_argument
00138     ,"sytrs(...) : Error, The number of rows in A and B must match."
00139     );
00140   FortranTypes::f_int info;
00141   LAPACK_Cpp::sytrs(  A.uplo(), A.rows(), B->cols(), A.gms().col_ptr(1)
00142     , A.gms().max_rows(), ipiv, B->col_ptr(1), B->max_rows()
00143     , &info );
00144   TEUCHOS_TEST_FOR_EXCEPTION(
00145     info < 0, std::invalid_argument
00146     ,"sytrs(...): Error, Invalid argument "
00147     << -info << " sent to LAPACK function xSYTRS(...)"
00148     );
00149 }
00150 
00151 void DenseLinAlgLAPack::getrf(
00152   DMatrixSlice* A, FortranTypes::f_int ipiv[], FortranTypes::f_int* rank
00153   )
00154 {
00155   FortranTypes::f_int info;
00156   LAPACK_Cpp::getrf(
00157     A->rows(), A->cols(), A->col_ptr(1), A->max_rows()
00158     ,ipiv, &info
00159     );
00160   *rank = my_min( A->rows(), A->cols() );
00161   TEUCHOS_TEST_FOR_EXCEPTION(
00162     info < 0, std::invalid_argument
00163     ,"getrf(...): Error, Invalid argument "
00164     << -info << " sent to LAPACK function xGETRF(...)" );
00165   if(info > 0)
00166     *rank = info - 1;
00167 }
00168 
00169 void DenseLinAlgLAPack::getrs(
00170   const DMatrixSlice& LU, const FortranTypes::f_int ipiv[], BLAS_Cpp::Transp transp
00171   ,DMatrixSlice* B
00172   )
00173 {
00174   TEUCHOS_TEST_FOR_EXCEPTION(
00175     (LU.rows() != LU.cols() || LU.rows() != B->rows() ), std::invalid_argument
00176     ,"getrs(...) : Error, A must be square and the number of rows in A and B must match."
00177     );
00178   FortranTypes::f_int info;
00179   LAPACK_Cpp::getrs(
00180     transp, LU.rows(), B->cols(), LU.col_ptr(1), LU.max_rows(), ipiv
00181     ,B->col_ptr(1), B->max_rows(), &info
00182     );
00183   TEUCHOS_TEST_FOR_EXCEPTION(
00184     info < 0, std::invalid_argument
00185     ,"getrs(...): Error, Invalid argument "
00186     << -info << " sent to LAPACK function xGETRS(...)" );
00187   // If we get here B is the solution to the linear system.
00188 }
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Defines