ROL
LinearAlgebra.hpp
Go to the documentation of this file.
1 // @HEADER
2 // *****************************************************************************
3 // Rapid Optimization Library (ROL) Package
4 //
5 // Copyright 2014 NTESS and the ROL contributors.
6 // SPDX-License-Identifier: BSD-3-Clause
7 // *****************************************************************************
8 // @HEADER
9 
10 #include"Teuchos_LAPACK.hpp"
11 #include"ROL_StdVector.hpp"
12 
13 #ifndef __LINEAR_ALGEBRA__
14 #define __LINEAR_ALGEBRA__
15 
23 template<class Real>
24 void trisolve(ROL::Ptr<Teuchos::LAPACK<int,Real> > lapack,
25  const std::vector<Real>& a,
26  const std::vector<Real>& b,
27  const std::vector<Real>& c,
28  const std::vector<Real>& r,
29  std::vector<Real>& x ) {
30 
31  const char TRANS = 'N';
32  const int N = b.size();
33  const int NRHS = 1;
34  int info;
35 
36  std::vector<Real> dl(a);
37  std::vector<Real> d(b);
38  std::vector<Real> du(c);
39  std::vector<Real> du2(N-2,0.0);
40  std::vector<int> ipiv(N);
41 
42  // Do matrix factorization, overwriting the LU bands
43  lapack->GTTRF(N,&dl[0],&d[0],&du[0],&du2[0],&ipiv[0],&info);
44 
45  x = r;
46 
47  // Solve the system using the banded LU factors
48  lapack->GTTRS(TRANS,N,NRHS,&dl[0],&d[0],&du[0],&du2[0],&ipiv[0],&x[0],N,&info);
49 }
50 
51 
58 template<class Real>
59 void lusolve(ROL::Ptr<Teuchos::LAPACK<int,Real> > lapack,
60  const std::vector<Real> &A,
61  const std::vector<Real> &B,
62  std::vector<Real> &X) {
63 
64  const char TRANS = 'N';
65  const int N2 = A.size();
66  const int N = sqrt(N2);
67  const int M = int(B.size()/N);
68  const int LDA = N;
69  int info;
70  std::vector<int> ipiv(N);
71  std::vector<Real> PLU(A);
72  X = B;
73 
74  // Do matrix factorization
75  lapack->GETRF(N,N,&PLU[0],LDA,&ipiv[0],&info);
76 
77  // Solve LU-factored system
78  lapack->GETRS(TRANS,N,M,&PLU[0],LDA,&ipiv[0],&X[0],LDA,&info);
79 }
80 
81 
82 
89 template<class Real>
90 void cholsolve(ROL::Ptr<Teuchos::LAPACK<int,Real> > lapack,
91  const std::vector<Real> &A,
92  const std::vector<Real> &B,
93  std::vector<Real> &X) {
94 
95  const char UPLO = 'U';
96  const int N2 = A.size();
97  const int N = sqrt(N2);
98  const int M = int(B.size()/N);
99  const int LDA = N;
100  int info;
101  std::vector<Real> C(A);
102  X = B;
103 
104  // Do Cholesky matrix factorization
105  lapack->POTRF(UPLO,N,&C[0],LDA,&info);
106 
107  // Solve Cholesky-factored system
108  lapack->POTRS(UPLO,N,&C[0],LDA,&X[0],LDA,&info);
109 }
110 
111 #endif
112 
113 
void trisolve(ROL::Ptr< Teuchos::LAPACK< int, Real > > lapack, const std::vector< Real > &a, const std::vector< Real > &b, const std::vector< Real > &c, const std::vector< Real > &r, std::vector< Real > &x)
Solve a tridiagonal system.
void cholsolve(ROL::Ptr< Teuchos::LAPACK< int, Real > > lapack, const std::vector< Real > &A, const std::vector< Real > &B, std::vector< Real > &X)
void lusolve(ROL::Ptr< Teuchos::LAPACK< int, Real > > lapack, const std::vector< Real > &A, const std::vector< Real > &B, std::vector< Real > &X)