source: ThirdParty/vmg/src/solver/dsysv.hpp@ 33a694

Action_Thermostats Add_AtomRandomPerturbation Add_RotateAroundBondAction Add_SelectAtomByNameAction Adding_Graph_to_ChangeBondActions Adding_MD_integration_tests Adding_StructOpt_integration_tests AutomationFragmentation_failures Candidate_v1.6.1 ChangeBugEmailaddress ChangingTestPorts ChemicalSpaceEvaluator Docu_Python_wait EmpiricalPotential_contain_HomologyGraph_documentation Enhance_userguide Enhanced_StructuralOptimization Enhanced_StructuralOptimization_continued Example_ManyWaysToTranslateAtom Exclude_Hydrogens_annealWithBondGraph Fix_ChronosMutex Fix_StatusMsg Fix_StepWorldTime_single_argument Fix_Verbose_Codepatterns ForceAnnealing_goodresults ForceAnnealing_oldresults ForceAnnealing_tocheck ForceAnnealing_with_BondGraph ForceAnnealing_with_BondGraph_continued ForceAnnealing_with_BondGraph_continued_betteresults ForceAnnealing_with_BondGraph_contraction-expansion GeometryObjects Gui_displays_atomic_force_velocity IndependentFragmentGrids_IntegrationTest JobMarket_RobustOnKillsSegFaults JobMarket_StableWorkerPool PythonUI_with_named_parameters QtGui_reactivate_TimeChanged_changes Recreated_GuiChecks RotateToPrincipalAxisSystem_UndoRedo StoppableMakroAction TremoloParser_IncreasedPrecision TremoloParser_MultipleTimesteps Ubuntu_1604_changes
Last change on this file since 33a694 was 7faa5c, checked in by Frederik Heber <heber@…>, 8 years ago

Merge commit 'de061d9d851257a04e924d4472df4523d33bb08b' as 'ThirdParty/vmg'

  • Property mode set to 100644
File size: 5.4 KB
Line 
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 dsysv.hpp
21 * @author Julian Iseringhausen <isering@ins.uni-bonn.de>
22 * @date Mon Apr 18 13:09:54 2011
23 *
24 * @brief Lapack solver for symmetric matrices
25 *
26 */
27
28#ifndef DSYSV_HPP_
29#define DSYSV_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#include "base/defs.hpp"
39
40namespace VMG
41{
42
43extern "C"
44{
45
46void FC_FUNC(dsysv, DSYSV) (const char* UPLO, const int* N, const int* NRHS,
47 double* A, const int* LDA, int* IPIV, double* B,
48 const int* LDB, double* WORK, const int* LWORK,
49 int* INFO);
50
51void FC_FUNC(dsyrfs, DSYRFS) (const char* UPLO, const int* N, const int* NRHS,
52 const vmg_float* A, const int* LDA, double* AF,
53 const int* LDAF, const int* IPIV, const vmg_float* B,
54 const int* LDB, double* X, const int* LDX,
55 double* FERR, double* BERR, double* WORK, int* IWORK,
56 int* INFO);
57
58void FC_FUNC(ssysv, SSYSV) (const char* UPLO, const int* N, const int* NRHS, float* A,
59 const int* LDA, int* IPIV, float* B, const int* LDB,
60 float* WORK, const int* LWORK, int* INFO);
61
62void FC_FUNC(ssyrfs, SSYRFS) (const char* UPLO, const int* N, const int* NRHS,
63 const vmg_float* A, const int* LDA, float* AF,
64 const int* LDAF, const int* IPIV, const vmg_float* B,
65 const int* LDB, float* X, const int* LDX, float* FERR,
66 float* BERR, double* WORK, int* IWORK, int* INFO);
67
68} /* extern "C" */
69
70/* By default, assume fcs_float is double. */
71#if defined(FCS_FLOAT_IS_FLOAT)
72#define dsysv FC_FUNC(ssysv, SSYSV)
73#define dsyrfs FC_FUNC(ssyrfs, SSYRFS)
74#else
75#define dsysv FC_FUNC(dsysv, DSYSV)
76#define dsyrfs FC_FUNC(dsyrfs, DSYRFS)
77#endif
78
79template<class T>
80class DSYSV : public T
81{
82public:
83 DSYSV(bool register_ = true) :
84 T(register_)
85 {Init();}
86
87 DSYSV(int size, bool register_ = true) :
88 T(size, register_)
89 {Init();}
90
91 virtual ~DSYSV();
92
93protected:
94 void Compute();
95
96private:
97 void Init();
98 void Realloc();
99
100 char la_uplo;
101 int la_n, la_nrhs, *la_ipiv, la_lwork, la_info, *la_iwork;
102 vmg_float *la_A, *la_A_orig, *la_B, *la_B_orig, *la_work, *la_work2;
103
104 int cur_lapack_size, max_lapack_size;
105};
106
107template<class T>
108void DSYSV<T>::Compute()
109{
110 vmg_float ferr, berr;
111 vmg_float opt_work;
112
113 this->Realloc();
114
115 for (int i=0; i<this->cur_lapack_size; i++) {
116 la_B[i] = la_B_orig[i] = this->Rhs(i);
117 for (int j=i; j<this->cur_lapack_size; j++)
118 la_A[j + i*this->cur_lapack_size] = la_A_orig[j + i*this->cur_lapack_size] = this->Mat(i,j);
119 }
120
121 // Determine optimal size of working space
122 this->la_lwork = -1;
123 dsysv (&la_uplo, &la_n, &la_nrhs, la_A, &la_n, la_ipiv, la_B, &la_n, &opt_work, &la_lwork, &la_info);
124
125 // Resize working space
126 this->la_lwork = static_cast<int>(opt_work);
127 delete [] this->la_work;
128 this->la_work = new vmg_float[this->la_lwork];
129
130 // Solve system
131 dsysv (&la_uplo, &la_n, &la_nrhs, la_A, &la_n, la_ipiv, la_B, &la_n, la_work, &la_lwork, &la_info);
132
133 // Improve computed solution
134 dsyrfs (&la_uplo, &la_n, &la_nrhs, la_A_orig, &la_n, la_A, &la_n, la_ipiv, la_B_orig, &la_n, la_B, &la_n, &ferr, &berr, la_work2, la_iwork, &la_info);
135
136 // Write solution back
137 for (int i=0; i<this->cur_lapack_size; i++)
138 this->Sol(i) = this->la_B[i];
139}
140
141template<class T>
142void DSYSV<T>::Realloc()
143{
144 this->cur_lapack_size = this->Size();
145
146 if (this->cur_lapack_size > this->max_lapack_size) {
147
148 delete [] la_A;
149 delete [] la_A_orig;
150 delete [] la_B;
151 delete [] la_B_orig;
152 delete [] la_ipiv;
153 delete [] la_work2;
154 delete [] la_iwork;
155
156 this->la_A = new vmg_float[this->cur_lapack_size * this->cur_lapack_size];
157 this->la_A_orig = new vmg_float[this->cur_lapack_size * this->cur_lapack_size];
158 this->la_B = new vmg_float[this->cur_lapack_size];
159 this->la_B_orig = new vmg_float[this->cur_lapack_size];
160 this->la_ipiv = new int[this->cur_lapack_size];
161 this->la_work2 = new vmg_float[3 * this->cur_lapack_size];
162 this->la_iwork = new int[this->cur_lapack_size];
163
164 this->la_n = this->cur_lapack_size;
165
166 this->max_lapack_size = this->cur_lapack_size;
167
168 }
169}
170
171template<class T>
172void DSYSV<T>::Init()
173{
174 this->cur_lapack_size = 0;
175 this->max_lapack_size = 0;
176
177 this->la_A = NULL;
178 this->la_A_orig = NULL;
179 this->la_B = NULL;
180 this->la_B_orig = NULL;
181 this->la_work = NULL;
182 this->la_ipiv = NULL;
183 this->la_work2 = NULL;
184 this->la_iwork = NULL;
185
186 this->la_nrhs = 1;
187 this->la_uplo = 'L';
188}
189
190template<class T>
191DSYSV<T>::~DSYSV()
192{
193 delete [] this->la_A;
194 delete [] this->la_A_orig;
195 delete [] this->la_B;
196 delete [] this->la_B_orig;
197 delete [] this->la_work;
198 delete [] this->la_ipiv;
199 delete [] this->la_work2;
200 delete [] this->la_iwork;
201}
202
203}
204
205#endif /* DSYSV_HPP_ */
Note: See TracBrowser for help on using the repository browser.