1 | /*
|
---|
2 | * vmg - a versatile multigrid solver
|
---|
3 | * Copyright (C) 2012 Institute for Numerical Simulation, University of Bonn
|
---|
4 | *
|
---|
5 | * vmg is free software: you can redistribute it and/or modify
|
---|
6 | * it under the terms of the GNU General Public License as published by
|
---|
7 | * the Free Software Foundation, either version 3 of the License, or
|
---|
8 | * (at your option) any later version.
|
---|
9 | *
|
---|
10 | * vmg is distributed in the hope that it will be useful,
|
---|
11 | * but WITHOUT ANY WARRANTY; without even the implied warranty of
|
---|
12 | * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
|
---|
13 | * GNU General Public License for more details.
|
---|
14 | *
|
---|
15 | * You should have received a copy of the GNU General Public License
|
---|
16 | * along with this program. If not, see <http://www.gnu.org/licenses/>.
|
---|
17 | */
|
---|
18 |
|
---|
19 | /**
|
---|
20 | * @file dgesv.hpp
|
---|
21 | * @author Julian Iseringhausen <isering@ins.uni-bonn.de>
|
---|
22 | * @date Mon Apr 18 13:09:29 2011
|
---|
23 | *
|
---|
24 | * @brief Lapack solver for general matrices.
|
---|
25 | *
|
---|
26 | */
|
---|
27 |
|
---|
28 | #ifndef DGESV_HPP_
|
---|
29 | #define DGESV_HPP_
|
---|
30 |
|
---|
31 | #ifndef HAVE_LAPACK
|
---|
32 | #error You need LAPACK in order to use MGDSYSV
|
---|
33 | #endif
|
---|
34 |
|
---|
35 | #include <cassert>
|
---|
36 | #include <cstddef>
|
---|
37 |
|
---|
38 | namespace VMG
|
---|
39 | {
|
---|
40 |
|
---|
41 | extern "C" {
|
---|
42 |
|
---|
43 | void FC_FUNC(dgesv,DGESV) (const int* N, const int* NRHS, vmg_float* A,
|
---|
44 | const int* LDA, int* IPIV, vmg_float* B,
|
---|
45 | const int* LDB, int* INFO);
|
---|
46 |
|
---|
47 | void FC_FUNC(sgesv, SGESV) (const int* N, const int* NRHS, vmg_float* A,
|
---|
48 | const int* LDA, int* IPIV, vmg_float* B,
|
---|
49 | const int* LDB, int* INFO);
|
---|
50 |
|
---|
51 | } /* extern "C" */
|
---|
52 |
|
---|
53 | /* By default, assume fcs_float is double. */
|
---|
54 | #if defined(FCS_FLOAT_IS_FLOAT)
|
---|
55 | #define dgesv FC_FUNC(sgesv,SGESV)
|
---|
56 | #else
|
---|
57 | #define dgesv FC_FUNC(dgesv,DGESV)
|
---|
58 | #endif
|
---|
59 |
|
---|
60 | template<class T>
|
---|
61 | class DGESV : public T
|
---|
62 | {
|
---|
63 | public:
|
---|
64 | DGESV(bool register_ = true) :
|
---|
65 | T(register_)
|
---|
66 | {Init();}
|
---|
67 |
|
---|
68 | DGESV(int size, bool register_ = true) :
|
---|
69 | T(size, register_)
|
---|
70 | {Init();}
|
---|
71 |
|
---|
72 | virtual ~DGESV();
|
---|
73 |
|
---|
74 | protected:
|
---|
75 | void Compute();
|
---|
76 |
|
---|
77 | private:
|
---|
78 | void Init();
|
---|
79 | void Realloc();
|
---|
80 |
|
---|
81 | int la_INFO, la_LDA, la_LDB, la_N, la_NRHS;
|
---|
82 | int *la_IPIV;
|
---|
83 | vmg_float *la_A, *la_B;
|
---|
84 |
|
---|
85 | int la_cur_size, la_max_size;
|
---|
86 | };
|
---|
87 |
|
---|
88 | template<class T>
|
---|
89 | void DGESV<T>::Compute()
|
---|
90 | {
|
---|
91 | // Adjust size of vectors
|
---|
92 | this->Realloc();
|
---|
93 |
|
---|
94 | la_N = la_LDA = la_LDB = la_cur_size;
|
---|
95 |
|
---|
96 | // Rewrite matrix in column-major order
|
---|
97 | for (int i=0; i<la_N; i++) {
|
---|
98 | la_B[i] = this->Rhs(i);
|
---|
99 | for (int j=0; j<la_N; j++)
|
---|
100 | la_A[i + j*la_N] = this->Mat(i,j);
|
---|
101 | }
|
---|
102 |
|
---|
103 | // Solve system of equations
|
---|
104 | dgesv(&la_N, &la_NRHS, la_A, &la_LDA, la_IPIV, la_B, &la_LDB, &la_INFO);
|
---|
105 |
|
---|
106 | // Assert successful exit of solver routine
|
---|
107 | assert(la_INFO == 0);
|
---|
108 |
|
---|
109 | // Write solution back
|
---|
110 | for (int i=0; i<la_N; i++)
|
---|
111 | this->Sol(i) = la_B[i];
|
---|
112 | }
|
---|
113 |
|
---|
114 | template<class T>
|
---|
115 | void DGESV<T>::Init()
|
---|
116 | {
|
---|
117 | la_cur_size = 0;
|
---|
118 | la_max_size = 0;
|
---|
119 | la_NRHS = 1;
|
---|
120 |
|
---|
121 | la_IPIV = NULL;
|
---|
122 | la_A = NULL;
|
---|
123 | la_B = NULL;
|
---|
124 | }
|
---|
125 |
|
---|
126 | template<class T>
|
---|
127 | void DGESV<T>::Realloc()
|
---|
128 | {
|
---|
129 | la_cur_size = this->Size();
|
---|
130 |
|
---|
131 | if (la_cur_size > la_max_size) {
|
---|
132 |
|
---|
133 | delete [] la_IPIV;
|
---|
134 | delete [] la_A;
|
---|
135 | delete [] la_B;
|
---|
136 |
|
---|
137 | la_IPIV = new int[la_cur_size];
|
---|
138 | la_A = new vmg_float[la_cur_size*la_cur_size];
|
---|
139 | la_B = new vmg_float[la_cur_size];
|
---|
140 |
|
---|
141 | la_max_size = la_cur_size;
|
---|
142 |
|
---|
143 | }
|
---|
144 | }
|
---|
145 |
|
---|
146 | template<class T>
|
---|
147 | DGESV<T>::~DGESV()
|
---|
148 | {
|
---|
149 | delete [] la_IPIV;
|
---|
150 | delete [] la_A;
|
---|
151 | delete [] la_B;
|
---|
152 | }
|
---|
153 |
|
---|
154 | }
|
---|
155 |
|
---|
156 | #endif /* DGESV_HPP_ */
|
---|