ROL
LinearAlgebra.hpp
Go to the documentation of this file.
1#include"Teuchos_LAPACK.hpp"
2#include"ROL_StdVector.hpp"
3
4#ifndef __LINEAR_ALGEBRA__
5#define __LINEAR_ALGEBRA__
6
14template<class Real>
15void trisolve(ROL::Ptr<Teuchos::LAPACK<int,Real> > lapack,
16 const std::vector<Real>& a,
17 const std::vector<Real>& b,
18 const std::vector<Real>& c,
19 const std::vector<Real>& r,
20 std::vector<Real>& x ) {
21
22 const char TRANS = 'N';
23 const int N = b.size();
24 const int NRHS = 1;
25 int info;
26
27 std::vector<Real> dl(a);
28 std::vector<Real> d(b);
29 std::vector<Real> du(c);
30 std::vector<Real> du2(N-2,0.0);
31 std::vector<int> ipiv(N);
32
33 // Do matrix factorization, overwriting the LU bands
34 lapack->GTTRF(N,&dl[0],&d[0],&du[0],&du2[0],&ipiv[0],&info);
35
36 x = r;
37
38 // Solve the system using the banded LU factors
39 lapack->GTTRS(TRANS,N,NRHS,&dl[0],&d[0],&du[0],&du2[0],&ipiv[0],&x[0],N,&info);
40}
41
42
49template<class Real>
50void lusolve(ROL::Ptr<Teuchos::LAPACK<int,Real> > lapack,
51 const std::vector<Real> &A,
52 const std::vector<Real> &B,
53 std::vector<Real> &X) {
54
55 const char TRANS = 'N';
56 const int N2 = A.size();
57 const int N = sqrt(N2);
58 const int M = int(B.size()/N);
59 const int LDA = N;
60 int info;
61 std::vector<int> ipiv(N);
62 std::vector<Real> PLU(A);
63 X = B;
64
65 // Do matrix factorization
66 lapack->GETRF(N,N,&PLU[0],LDA,&ipiv[0],&info);
67
68 // Solve LU-factored system
69 lapack->GETRS(TRANS,N,M,&PLU[0],LDA,&ipiv[0],&X[0],LDA,&info);
70}
71
72
73
80template<class Real>
81void cholsolve(ROL::Ptr<Teuchos::LAPACK<int,Real> > lapack,
82 const std::vector<Real> &A,
83 const std::vector<Real> &B,
84 std::vector<Real> &X) {
85
86 const char UPLO = 'U';
87 const int N2 = A.size();
88 const int N = sqrt(N2);
89 const int M = int(B.size()/N);
90 const int LDA = N;
91 int info;
92 std::vector<Real> C(A);
93 X = B;
94
95 // Do Cholesky matrix factorization
96 lapack->POTRF(UPLO,N,&C[0],LDA,&info);
97
98 // Solve Cholesky-factored system
99 lapack->POTRS(UPLO,N,&C[0],LDA,&X[0],LDA,&info);
100}
101
102#endif
103
104
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)
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.