Revisión | dbc6d52c4017fa9a06f2e5a16dec29710294530a (tree) |
---|---|
Tiempo | 2012-12-29 14:25:10 |
Autor | Katsuhiko Nishimra <ktns.87@gmai...> |
Commiter | Katsuhiko Nishimra |
Remove Blas/Lapack_GNU.cpp. #30409
git-svn-id: https://svn.sourceforge.jp/svnroot/molds/trunk@1221 1136aad2-a195-0410-b898-f5ea1d11b9d8
@@ -1,296 +0,0 @@ | ||
1 | -//************************************************************************// | |
2 | -// Copyright (C) 2011-2012 Mikiya Fujii // | |
3 | -// Copyright (C) 2012-2012 Katsuhiko Nishimra // | |
4 | -// // | |
5 | -// This file is part of MolDS. // | |
6 | -// // | |
7 | -// MolDS is free software: you can redistribute it and/or modify // | |
8 | -// it under the terms of the GNU General Public License as published by // | |
9 | -// the Free Software Foundation, either version 3 of the License, or // | |
10 | -// (at your option) any later version. // | |
11 | -// // | |
12 | -// MolDS is distributed in the hope that it will be useful, // | |
13 | -// but WITHOUT ANY WARRANTY; without even the implied warranty of // | |
14 | -// MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the // | |
15 | -// GNU General Public License for more details. // | |
16 | -// // | |
17 | -// You should have received a copy of the GNU General Public License // | |
18 | -// along with MolDS. If not, see <http://www.gnu.org/licenses/>. // | |
19 | -//************************************************************************// | |
20 | -#include<stdio.h> | |
21 | -#include<stdlib.h> | |
22 | -#include<iostream> | |
23 | -#include<sstream> | |
24 | -#include<math.h> | |
25 | -#include<string> | |
26 | -#include<stdexcept> | |
27 | -#include<boost/format.hpp> | |
28 | -#include"cblas.h" | |
29 | -#include"../base/PrintController.h" | |
30 | -#include"../base/MolDSException.h" | |
31 | -#include"../base/Uncopyable.h" | |
32 | -#include"Blas.h" | |
33 | -using namespace std; | |
34 | -using namespace MolDS_base; | |
35 | - | |
36 | -namespace MolDS_wrappers{ | |
37 | -Blas* Blas::blas = NULL; | |
38 | - | |
39 | -Blas::Blas(){ | |
40 | -} | |
41 | - | |
42 | -Blas::~Blas(){ | |
43 | -} | |
44 | - | |
45 | -Blas* Blas::GetInstance(){ | |
46 | - if(blas == NULL){ | |
47 | - blas = new Blas(); | |
48 | - //this->OutputLog("Blas created.\n\n"); | |
49 | - } | |
50 | - return blas; | |
51 | -} | |
52 | - | |
53 | -void Blas::DeleteInstance(){ | |
54 | - if(blas != NULL){ | |
55 | - delete blas; | |
56 | - //this->OutputLog("Blas deleted\n\n"); | |
57 | - } | |
58 | - blas = NULL; | |
59 | -} | |
60 | - | |
61 | -// vectorY = vectorX | |
62 | -// vectorX: n-vector | |
63 | -// vectorY: n-vector | |
64 | -void Blas::Dcopy(molds_blas_int n, | |
65 | - double const* vectorX, | |
66 | - double * vectorY)const{ | |
67 | - molds_blas_int incrementX=1; | |
68 | - molds_blas_int incrementY=1; | |
69 | - this->Dcopy(n, vectorX, incrementX, vectorY, incrementY); | |
70 | -} | |
71 | - | |
72 | -// vectorY = vectorX | |
73 | -// vectorX: n-vector | |
74 | -// vectorY: n-vector | |
75 | -void Blas::Dcopy(molds_blas_int n, | |
76 | - double const* vectorX, molds_blas_int incrementX, | |
77 | - double* vectorY, molds_blas_int incrementY) const{ | |
78 | - double* x = const_cast<double*>(&vectorX[0]); | |
79 | - cblas_dcopy(n, x, incrementX, vectorY, incrementY); | |
80 | -} | |
81 | - | |
82 | -// vectorY = alpha*vectorX + vectorY | |
83 | -// vectorX: n-vector | |
84 | -// vectorY: n-vector | |
85 | -void Blas::Daxpy(molds_blas_int n, double alpha, | |
86 | - double const* vectorX, | |
87 | - double* vectorY) const{ | |
88 | - molds_blas_int incrementX=1; | |
89 | - molds_blas_int incrementY=1; | |
90 | - this->Daxpy(n, alpha, vectorX, incrementX, vectorY, incrementY); | |
91 | -} | |
92 | - | |
93 | -// vectorY = alpha*vectorX + vectorY | |
94 | -// vectorX: n-vector | |
95 | -// vectorY: n-vector | |
96 | -void Blas::Daxpy(molds_blas_int n, double alpha, | |
97 | - double const* vectorX, molds_blas_int incrementX, | |
98 | - double* vectorY, molds_blas_int incrementY) const{ | |
99 | - double* x = const_cast<double*>(&vectorX[0]); | |
100 | - cblas_daxpy(n, alpha, x, incrementX, vectorY, incrementY); | |
101 | -} | |
102 | - | |
103 | -// returns vectorX^T*vectorY | |
104 | -// vectorX: n-vector | |
105 | -// vectorY: n-vector | |
106 | -double Blas::Ddot(molds_blas_int n, | |
107 | - double const* vectorX, | |
108 | - double const* vectorY) const{ | |
109 | - molds_blas_int incrementX=1; | |
110 | - molds_blas_int incrementY=1; | |
111 | - return this->Ddot(n, vectorX, incrementX, vectorY, incrementY); | |
112 | -} | |
113 | - | |
114 | -// returns vectorX^T*vectorY | |
115 | -// vectorX: n-vector | |
116 | -// vectorY: n-vector | |
117 | -double Blas::Ddot(molds_blas_int n, | |
118 | - double const* vectorX, molds_blas_int incrementX, | |
119 | - double const* vectorY, molds_blas_int incrementY)const{ | |
120 | - double* x=const_cast<double*>(vectorX), | |
121 | - * y=const_cast<double*>(vectorY); | |
122 | - return cblas_ddot(n, x, incrementX, y, incrementY); | |
123 | -} | |
124 | - | |
125 | -// vectorY = matrixA*vectorX | |
126 | -// matrixA: m*n-matrix (matrixA[m][n] in row-major (C/C++ style)) | |
127 | -// vectorX: n-vector | |
128 | -// vectorY: m-vector | |
129 | -void Blas::Dgemv(molds_blas_int m, molds_blas_int n, | |
130 | - double const* const* matrixA, | |
131 | - double const* vectorX, | |
132 | - double* vectorY) const{ | |
133 | - bool isColumnMajorMatrixA = false; // because, in general, C/C++ style is row-major. | |
134 | - molds_blas_int incrementX=1; | |
135 | - molds_blas_int incrementY=1; | |
136 | - double alpha =1.0; | |
137 | - double beta =0.0; | |
138 | - this->Dgemv(isColumnMajorMatrixA, m, n, alpha, matrixA, vectorX, incrementX, beta, vectorY, incrementY); | |
139 | -} | |
140 | - | |
141 | -// vectorY = alpha*matrixA*vectorX + beta*vectorY | |
142 | -// matrixA: m*n-matrix | |
143 | -// vectorX: n-vector | |
144 | -// vectorY: m-vector | |
145 | -void Blas::Dgemv(bool isColumnMajorMatrixA, | |
146 | - molds_blas_int m, molds_blas_int n, | |
147 | - double alpha, | |
148 | - double const* const* matrixA, | |
149 | - double const* vectorX, | |
150 | - molds_blas_int incrementX, | |
151 | - double beta, | |
152 | - double* vectorY, | |
153 | - molds_blas_int incrementY) const{ | |
154 | - double* a = const_cast<double*>(&matrixA[0][0]); | |
155 | - double* x = const_cast<double*>(&vectorX[0]); | |
156 | - CBLAS_TRANSPOSE transA; | |
157 | - if(isColumnMajorMatrixA){ | |
158 | - transA = CblasNoTrans; | |
159 | - } | |
160 | - else{ | |
161 | - transA = CblasTrans; | |
162 | - swap(m,n); | |
163 | - } | |
164 | - molds_blas_int lda = m; | |
165 | - cblas_dgemv(CblasColMajor, transA, m, n, alpha, a, lda, x, incrementX, beta, vectorY, incrementY); | |
166 | - | |
167 | -} | |
168 | - | |
169 | -// vectorY = matrixA*vectorX | |
170 | -// matrixA: n*n-matrix,symmetric (Use the upper triangular part) | |
171 | -// vectorX: n-vector | |
172 | -// vectorY: n-vector | |
173 | -void Blas::Dsymv(molds_blas_int n, | |
174 | - double const* const* matrixA, | |
175 | - double const* vectorX, | |
176 | - double* vectorY) const{ | |
177 | - bool isColumnMajorMatrixA = false; // because, in general, C/C++ style is row-major. | |
178 | - molds_blas_int incrementX=1, incrementY=1; | |
179 | - double alpha=1.0, beta=0.0; | |
180 | - this->Dsymv(n, alpha, matrixA, vectorX, incrementX, beta, vectorY, incrementY); | |
181 | -} | |
182 | - | |
183 | -// vectorY = alpha*matrixA*vectorX + beta*vectorY | |
184 | -// matrixA: n*n-matrix,symmetric (Use the upper triangular part) | |
185 | -// vectorX: n-vector | |
186 | -// vectorY: n-vector | |
187 | -void Blas::Dsymv(molds_blas_int n, double alpha, | |
188 | - double const* const* matrixA, | |
189 | - double const* vectorX, molds_blas_int incrementX, | |
190 | - double beta, | |
191 | - double* vectorY, molds_blas_int incrementY) const{ | |
192 | - double* a = const_cast<double*>(&matrixA[0][0]); | |
193 | - double* x = const_cast<double*>(&vectorX[0]); | |
194 | - CBLAS_UPLO uploA=CblasUpper; | |
195 | - molds_blas_int lda = n; | |
196 | - cblas_dsymv(CblasRowMajor, uploA, n, alpha, a, lda, x, incrementX, beta, vectorY, incrementY); | |
197 | -} | |
198 | - | |
199 | -// matrixA = alpha*vectorX*vectorX^T + matrixA | |
200 | -// matrixA: n*n-matrix,symmetric (Use the upper triangular part, and copy it to the lower part.) | |
201 | -// vectorX: n-matrix | |
202 | -void Blas::Dsyr(molds_blas_int n, double alpha, | |
203 | - double const* vectorX, | |
204 | - double ** matrixA)const{ | |
205 | - molds_blas_int incrementX=1; | |
206 | - this->Dsyr(n, alpha, vectorX, incrementX, matrixA); | |
207 | -} | |
208 | - | |
209 | -void Blas::Dsyr(molds_blas_int n, double alpha, | |
210 | - double const* vectorX, molds_blas_int incrementX, | |
211 | - double ** matrixA)const{ | |
212 | - double* a = &matrixA[0][0]; | |
213 | - double* x = const_cast<double*>(&vectorX[0]); | |
214 | - CBLAS_UPLO uploA=CblasUpper; | |
215 | - molds_blas_int lda = n; | |
216 | - cblas_dsyr(CblasRowMajor, uploA, n, alpha, x, incrementX, a, lda); | |
217 | -#pragma omp parallel for schedule(auto) | |
218 | - for(molds_blas_int i=0;i<n;i++){ | |
219 | - for(molds_blas_int j=i+1;j<n;j++){ | |
220 | - matrixA[j][i] = matrixA[i][j]; | |
221 | - } | |
222 | - } | |
223 | -} | |
224 | - | |
225 | -// matrixC = matrixA*matrixB | |
226 | -// matrixA: m*k-matrix (matrixA[m][k] in row-major (C/C++ style)) | |
227 | -// matrixB: k*n-matrix (matrixB[k][n] in row-major (C/C++ style)) | |
228 | -// matrixC: m*n-matrix (matrixC[m][n] in row-major (C/C++ style)) | |
229 | -void Blas::Dgemm(molds_blas_int m, molds_blas_int n, molds_blas_int k, | |
230 | - double const* const* matrixA, | |
231 | - double const* const* matrixB, | |
232 | - double** matrixC) const{ | |
233 | - bool isColumnMajorMatrixA = false; // because, in general, C/C++ style is row-major. | |
234 | - bool isColumnMajorMatrixB = false; // because, in general, C/C++ style is row-major. | |
235 | - double alpha=1.0; | |
236 | - double beta =0.0; | |
237 | - this->Dgemm(isColumnMajorMatrixA, isColumnMajorMatrixB, m, n, k, alpha, matrixA, matrixB, beta, matrixC); | |
238 | -} | |
239 | - | |
240 | -// matrixC = alpha*matrixA*matrixB + beta*matrixC | |
241 | -// matrixA: m*k-matrix | |
242 | -// matrixB: k*n-matrix | |
243 | -// matrixC: m*n-matrix (matrixC[m][n] in row-major (C/C++ style)) | |
244 | -void Blas::Dgemm(bool isColumnMajorMatrixA, | |
245 | - bool isColumnMajorMatrixB, | |
246 | - molds_blas_int m, molds_blas_int n, molds_blas_int k, | |
247 | - double alpha, | |
248 | - double const* const* matrixA, | |
249 | - double const* const* matrixB, | |
250 | - double beta, | |
251 | - double** matrixC) const{ | |
252 | - double* a = const_cast<double*>(&matrixA[0][0]); | |
253 | - double* b = const_cast<double*>(&matrixB[0][0]); | |
254 | - double* c = &matrixC[0][0]; | |
255 | - | |
256 | - molds_blas_int lda; | |
257 | - CBLAS_TRANSPOSE transA; | |
258 | - if(isColumnMajorMatrixA){ | |
259 | - transA = CblasNoTrans; | |
260 | - lda = m; | |
261 | - } | |
262 | - else{ | |
263 | - transA = CblasTrans; | |
264 | - lda = k; | |
265 | - } | |
266 | - | |
267 | - molds_blas_int ldb; | |
268 | - CBLAS_TRANSPOSE transB; | |
269 | - if(isColumnMajorMatrixB){ | |
270 | - transB = CblasNoTrans; | |
271 | - ldb = k; | |
272 | - } | |
273 | - else{ | |
274 | - transB = CblasTrans; | |
275 | - ldb = n; | |
276 | - } | |
277 | - | |
278 | - double* tmpC; | |
279 | - tmpC = (double*)malloc( sizeof(double)*m*n); | |
280 | - for(molds_blas_int i=0; i<m; i++){ | |
281 | - for(molds_blas_int j=0; j<n; j++){ | |
282 | - tmpC[i+j*m] = matrixC[i][j]; | |
283 | - } | |
284 | - } | |
285 | - molds_blas_int ldc = m; | |
286 | - //call blas | |
287 | - cblas_dgemm(CblasColMajor, transA, transB, m, n, k, alpha, a, lda, b, ldb, beta, tmpC, ldc); | |
288 | - for(molds_blas_int i=0; i<m; i++){ | |
289 | - for(molds_blas_int j=0; j<n; j++){ | |
290 | - matrixC[i][j] = tmpC[i+j*m]; | |
291 | - } | |
292 | - } | |
293 | - free(tmpC); | |
294 | -} | |
295 | - | |
296 | -} |
@@ -1,365 +0,0 @@ | ||
1 | -//************************************************************************// | |
2 | -// Copyright (C) 2011-2012 Mikiya Fujii // | |
3 | -// // | |
4 | -// This file is part of MolDS. // | |
5 | -// // | |
6 | -// MolDS is free software: you can redistribute it and/or modify // | |
7 | -// it under the terms of the GNU General Public License as published by // | |
8 | -// the Free Software Foundation, either version 3 of the License, or // | |
9 | -// (at your option) any later version. // | |
10 | -// // | |
11 | -// MolDS is distributed in the hope that it will be useful, // | |
12 | -// but WITHOUT ANY WARRANTY; without even the implied warranty of // | |
13 | -// MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the // | |
14 | -// GNU General Public License for more details. // | |
15 | -// // | |
16 | -// You should have received a copy of the GNU General Public License // | |
17 | -// along with MolDS. If not, see <http://www.gnu.org/licenses/>. // | |
18 | -//************************************************************************// | |
19 | -#include<stdio.h> | |
20 | -#include<stdlib.h> | |
21 | -#include<limits.h> | |
22 | -#include<iostream> | |
23 | -#include<sstream> | |
24 | -#include<math.h> | |
25 | -#include<string> | |
26 | -#include<stdexcept> | |
27 | -#include<boost/format.hpp> | |
28 | -#if ( __WORDSIZE == 32 ) | |
29 | -#else | |
30 | -#define HAVE_LAPACK_CONFIG_H | |
31 | -#define LAPACK_ILP64 | |
32 | -#endif | |
33 | -#include"lapacke.h" | |
34 | -#include"../base/PrintController.h" | |
35 | -#include"../base/MolDSException.h" | |
36 | -#include"../base/Uncopyable.h" | |
37 | -#include"Lapack.h" | |
38 | -using namespace std; | |
39 | -using namespace MolDS_base; | |
40 | - | |
41 | -namespace MolDS_wrappers{ | |
42 | -Lapack* Lapack::lapack = NULL; | |
43 | - | |
44 | -Lapack::Lapack(){ | |
45 | - this->errorMessageDsyevdInfo = "Error in wrappers::Lapack::Dsyevd: info != 0: info = "; | |
46 | - this->errorMessageDsyevdSize = "Error in wrappers::Lapack::Dsyevd: size of matirx < 1\n"; | |
47 | - this->errorMessageDsysvInfo = "Error in wrappers::Lapack::Dsysv: info != 0: info = "; | |
48 | - this->errorMessageDsysvSize = "Error in wrappers::Lapack::Dsysv: size of matirx < 1\n"; | |
49 | - this->errorMessageDgetrsInfo = "Error in wrappers::Lapack::Dgetrs: info != 0: info = "; | |
50 | - this->errorMessageDgetrsSize = "Error in wrappers::Lapack::Dgetrs: size of matirx < 1\n"; | |
51 | - this->errorMessageDgetrfInfo = "Error in wrappers::Lapack::Dgetrf: info != 0: info = "; | |
52 | -} | |
53 | - | |
54 | -Lapack::~Lapack(){ | |
55 | -} | |
56 | - | |
57 | -Lapack* Lapack::GetInstance(){ | |
58 | - if(lapack == NULL){ | |
59 | - lapack = new Lapack(); | |
60 | - //this->OutputLog("Lapack created.\n\n"); | |
61 | - } | |
62 | - return lapack; | |
63 | -} | |
64 | - | |
65 | -void Lapack::DeleteInstance(){ | |
66 | - if(lapack != NULL){ | |
67 | - delete lapack; | |
68 | - //this->OutputLog("Lapack deleted\n\n"); | |
69 | - } | |
70 | - lapack = NULL; | |
71 | -} | |
72 | - | |
73 | - | |
74 | -/*** | |
75 | - * | |
76 | - * return notice | |
77 | - * i-th eigen value is eigenValues[i]. | |
78 | - * i-th eigen vector is (matirx[i][0], matirx[i][1], matirx[i][2], ....). | |
79 | - * | |
80 | - * ***/ | |
81 | -molds_lapack_int Lapack::Dsyevd(double** matrix, double* eigenValues, molds_lapack_int size, bool calcEigenVectors){ | |
82 | - molds_lapack_int info = 0; | |
83 | - molds_lapack_int k = 0; | |
84 | - molds_lapack_int lwork; | |
85 | - molds_lapack_int liwork; | |
86 | - char job; | |
87 | - char uplo = 'U'; | |
88 | - molds_lapack_int lda = size; | |
89 | - double* convertedMatrix; | |
90 | - double* tempEigenValues; | |
91 | - double* work; | |
92 | - molds_lapack_int* iwork; | |
93 | - | |
94 | - // set job type | |
95 | - if(calcEigenVectors){ | |
96 | - job = 'V'; | |
97 | - } | |
98 | - else{ | |
99 | - job = 'N'; | |
100 | - } | |
101 | - | |
102 | - // calc. lwork and liwork | |
103 | - if(size < 1 ){ | |
104 | - stringstream ss; | |
105 | - ss << errorMessageDsyevdSize; | |
106 | - throw MolDSException(ss.str()); | |
107 | - } | |
108 | - else if(size == 1){ | |
109 | - lwork = 1; | |
110 | - liwork = 1; | |
111 | - } | |
112 | - else if(1 < size && job == 'N'){ | |
113 | - lwork = 2*size + 1; | |
114 | - liwork = 2; | |
115 | - } | |
116 | - else{ | |
117 | - // calc. k | |
118 | - double temp = log((double)size)/log(2.0); | |
119 | - if( (double)((molds_lapack_int)temp) < temp ){ | |
120 | - k = (molds_lapack_int)temp + 1; | |
121 | - } | |
122 | - else{ | |
123 | - k = (molds_lapack_int)temp; | |
124 | - } | |
125 | - lwork = 3*size*size + (5+2*k)*size + 1; | |
126 | - liwork = 5*size + 3; | |
127 | - } | |
128 | - | |
129 | - // malloc | |
130 | - work = (double*)LAPACKE_malloc( sizeof(double)*lwork ); | |
131 | - iwork = (molds_lapack_int*)LAPACKE_malloc( sizeof(molds_lapack_int)*liwork ); | |
132 | - convertedMatrix = (double*)LAPACKE_malloc( sizeof(double)*size*size ); | |
133 | - tempEigenValues = (double*)LAPACKE_malloc( sizeof(double)*size ); | |
134 | - | |
135 | - for(molds_lapack_int i = 0; i < size; i++){ | |
136 | - for(molds_lapack_int j = i; j < size; j++){ | |
137 | - convertedMatrix[i+j*size] = matrix[i][j]; | |
138 | - } | |
139 | - } | |
140 | - | |
141 | - // call Lapack | |
142 | - info = LAPACKE_dsyevd_work(LAPACK_COL_MAJOR, job, uplo, size, convertedMatrix, lda, tempEigenValues, work, lwork, iwork, liwork); | |
143 | - | |
144 | - for(molds_lapack_int i = 0; i < size; i++){ | |
145 | - for(molds_lapack_int j = 0; j < size; j++){ | |
146 | - matrix[i][j] = convertedMatrix[j+i*size]; //i-th row is i-th eigen vector | |
147 | - //matrix[j][i] = convertedMatrix[j+i*size]; //i-th column is i-th eigen vector | |
148 | - } | |
149 | - } | |
150 | - | |
151 | - for(molds_lapack_int i=0;i<size;i++){ | |
152 | - double temp = 0.0; | |
153 | - for(molds_lapack_int j=0;j<size;j++){ | |
154 | - temp += matrix[i][j]; | |
155 | - } | |
156 | - if(temp<0){ | |
157 | - for(molds_lapack_int j=0;j<size;j++){ | |
158 | - matrix[i][j]*=-1.0; | |
159 | - } | |
160 | - } | |
161 | - } | |
162 | - | |
163 | - for(molds_lapack_int i = 0; i < size; i++){ | |
164 | - eigenValues[i] = tempEigenValues[i]; | |
165 | - } | |
166 | - //this->OutputLog(boost::format("size=%d lwork=%d liwork=%d k=%d info=%d\n") % size % lwork % liwork % k % info); | |
167 | - | |
168 | - // free | |
169 | - LAPACKE_free(work); | |
170 | - LAPACKE_free(iwork); | |
171 | - LAPACKE_free(convertedMatrix); | |
172 | - LAPACKE_free(tempEigenValues); | |
173 | - | |
174 | - if(info != 0){ | |
175 | - stringstream ss; | |
176 | - ss << errorMessageDsyevdInfo; | |
177 | - ss << info << endl; | |
178 | - throw MolDSException(ss.str()); | |
179 | - } | |
180 | - return info; | |
181 | -} | |
182 | - | |
183 | -/*** | |
184 | - * | |
185 | - * Slove matrix*X=b, then we get X by this method. | |
186 | - * The X is stored in b. | |
187 | - * | |
188 | - */ | |
189 | -molds_lapack_int Lapack::Dsysv(double const* const* matrix, double* b, molds_lapack_int size){ | |
190 | - molds_lapack_int info = 0; | |
191 | - molds_lapack_int lwork; | |
192 | - char uplo = 'U'; | |
193 | - molds_lapack_int lda = size; | |
194 | - molds_lapack_int ldb = size; | |
195 | - molds_lapack_int nrhs = 1; | |
196 | - double* convertedMatrix; | |
197 | - double* work; | |
198 | - double* tempB; | |
199 | - molds_lapack_int* ipiv; | |
200 | - | |
201 | - if(size < 1 ){ | |
202 | - stringstream ss; | |
203 | - ss << errorMessageDsysvSize; | |
204 | - throw MolDSException(ss.str()); | |
205 | - } | |
206 | - | |
207 | - // malloc | |
208 | - ipiv = (molds_lapack_int*)LAPACKE_malloc( sizeof(molds_lapack_int)*2*size ); | |
209 | - convertedMatrix = (double*)LAPACKE_malloc( sizeof(double)*size*size ); | |
210 | - tempB = (double*)LAPACKE_malloc( sizeof(double)*size ); | |
211 | - | |
212 | - for(molds_lapack_int i = 0; i < size; i++){ | |
213 | - for(molds_lapack_int j = i; j < size; j++){ | |
214 | - convertedMatrix[i+j*size] = matrix[i][j]; | |
215 | - } | |
216 | - } | |
217 | - for(molds_lapack_int i = 0; i < size; i++){ | |
218 | - tempB[i] = b[i]; | |
219 | - } | |
220 | - | |
221 | - // calc. lwork | |
222 | - double blockSize=0.0; | |
223 | -#pragma omp critical | |
224 | - { | |
225 | - lwork = -1; | |
226 | - double tempWork[3]={0.0, 0.0, 0.0}; | |
227 | - info = LAPACKE_dsysv_work(LAPACK_COL_MAJOR, uplo, size, nrhs, convertedMatrix, lda, ipiv, tempB, ldb, tempWork, lwork); | |
228 | - blockSize = tempWork[0]/size; | |
229 | - } | |
230 | - info = 0; | |
231 | - lwork = blockSize*size; | |
232 | - work = (double*)LAPACKE_malloc( sizeof(double)*lwork ); | |
233 | - | |
234 | - // call Lapack | |
235 | - info = LAPACKE_dsysv_work(LAPACK_COL_MAJOR, uplo, size, nrhs, convertedMatrix, lda, ipiv, tempB, ldb, work, lwork); | |
236 | - for(molds_lapack_int i = 0; i < size; i++){ | |
237 | - b[i] = tempB[i]; | |
238 | - } | |
239 | - | |
240 | - // free | |
241 | - LAPACKE_free(convertedMatrix); | |
242 | - LAPACKE_free(ipiv); | |
243 | - LAPACKE_free(work); | |
244 | - LAPACKE_free(tempB); | |
245 | - | |
246 | - if(info != 0){ | |
247 | - stringstream ss; | |
248 | - ss << errorMessageDsysvInfo; | |
249 | - ss << info << endl; | |
250 | - throw MolDSException(ss.str()); | |
251 | - } | |
252 | - return info; | |
253 | -} | |
254 | - | |
255 | -/*** | |
256 | - * | |
257 | - * Slove matrix*X[i]=b[i] (i=0, 1, ... , nrhs-1), then we get X[i] by this method. | |
258 | - * The X[i] is stored in b[i]. | |
259 | - * b[i][j] is j-th element of i-th solution, b[i]. | |
260 | - * | |
261 | - */ | |
262 | -molds_lapack_int Lapack::Dgetrs(double const* const* matrix, double** b, molds_lapack_int size, molds_lapack_int nrhs) const{ | |
263 | - molds_lapack_int info = 0; | |
264 | - char trans = 'N'; | |
265 | - molds_lapack_int lda = size; | |
266 | - molds_lapack_int ldb = size; | |
267 | - double* convertedMatrix; | |
268 | - double* convertedB; | |
269 | - molds_lapack_int* ipiv; | |
270 | - | |
271 | - if(size < 1 ){ | |
272 | - stringstream ss; | |
273 | - ss << errorMessageDgetrsSize; | |
274 | - throw MolDSException(ss.str()); | |
275 | - } | |
276 | - | |
277 | - try{ | |
278 | - // malloc | |
279 | - ipiv = (molds_lapack_int*)LAPACKE_malloc( sizeof(molds_lapack_int)*2*size); | |
280 | - convertedMatrix = (double*)LAPACKE_malloc( sizeof(double)*size*size); | |
281 | - convertedB = (double*)LAPACKE_malloc( sizeof(double)*nrhs*size); | |
282 | - for(molds_lapack_int i = 0; i < size; i++){ | |
283 | - for(molds_lapack_int j = 0; j < size; j++){ | |
284 | - convertedMatrix[i+j*size] = matrix[i][j]; | |
285 | - } | |
286 | - } | |
287 | - for(molds_lapack_int i = 0; i < nrhs; i++){ | |
288 | - for(molds_lapack_int j = 0; j < size; j++){ | |
289 | - convertedB[j+i*size] = b[i][j]; | |
290 | - } | |
291 | - } | |
292 | - this->Dgetrf(convertedMatrix, ipiv, size, size); | |
293 | - info = LAPACKE_dgetrs_work(LAPACK_COL_MAJOR, trans, size, nrhs, convertedMatrix, lda, ipiv, convertedB, ldb); | |
294 | - for(molds_lapack_int i = 0; i < nrhs; i++){ | |
295 | - for(molds_lapack_int j = 0; j < size; j++){ | |
296 | - b[i][j] = convertedB[j+i*size]; | |
297 | - } | |
298 | - } | |
299 | - } | |
300 | - catch(MolDSException ex){ | |
301 | - // free | |
302 | - LAPACKE_free(convertedMatrix); | |
303 | - LAPACKE_free(convertedB); | |
304 | - LAPACKE_free(ipiv); | |
305 | - throw ex; | |
306 | - } | |
307 | - // free | |
308 | - LAPACKE_free(convertedMatrix); | |
309 | - LAPACKE_free(convertedB); | |
310 | - LAPACKE_free(ipiv); | |
311 | - | |
312 | - if(info != 0){ | |
313 | - stringstream ss; | |
314 | - ss << errorMessageDgetrsInfo; | |
315 | - ss << info << endl; | |
316 | - throw MolDSException(ss.str()); | |
317 | - } | |
318 | - return info; | |
319 | -} | |
320 | - | |
321 | -// Argument "matrix" is sizeM*sizeN matrix. | |
322 | -// Argument "matrix" will be LU-decomposed. | |
323 | -molds_lapack_int Lapack::Dgetrf(double** matrix, molds_lapack_int sizeM, molds_lapack_int sizeN) const{ | |
324 | - molds_lapack_int* ipiv = (molds_lapack_int*)LAPACKE_malloc( sizeof(molds_lapack_int)*2*sizeM ); | |
325 | - this->Dgetrf(matrix, ipiv, sizeM, sizeN); | |
326 | - LAPACKE_free(ipiv); | |
327 | - molds_lapack_int info = 0; | |
328 | - return info; | |
329 | -} | |
330 | - | |
331 | -// Argument "matrix" is sizeM*sizeN matrix. | |
332 | -// Argument "matrix" will be LU-decomposed. | |
333 | -molds_lapack_int Lapack::Dgetrf(double** matrix, molds_lapack_int* ipiv, molds_lapack_int sizeM, molds_lapack_int sizeN) const{ | |
334 | - double* convertedMatrix = (double*)LAPACKE_malloc( sizeof(double)*sizeM*sizeN ); | |
335 | - for(molds_lapack_int i=0; i<sizeM; i++){ | |
336 | - for(molds_lapack_int j=0; j<sizeN; j++){ | |
337 | - convertedMatrix[i+j*sizeM] = matrix[i][j]; | |
338 | - } | |
339 | - } | |
340 | - this->Dgetrf(convertedMatrix, ipiv, sizeM, sizeN); | |
341 | - for(molds_lapack_int i=0; i<sizeM; i++){ | |
342 | - for(molds_lapack_int j=0; j<sizeN; j++){ | |
343 | - matrix[i][j] = convertedMatrix[i+j*sizeM]; | |
344 | - } | |
345 | - } | |
346 | - LAPACKE_free(convertedMatrix); | |
347 | - molds_lapack_int info = 0; | |
348 | - return info; | |
349 | -} | |
350 | - | |
351 | -// Argument "matrix" means sizeM * sizeN matrix. | |
352 | -// The each element of "matrix" should be stored in 1-dimensional vecotre with column major (Fortran type). | |
353 | -molds_lapack_int Lapack::Dgetrf(double* matrix, molds_lapack_int* ipiv, molds_lapack_int sizeM, molds_lapack_int sizeN) const{ | |
354 | - molds_lapack_int info = 0; | |
355 | - molds_lapack_int lda = sizeM; | |
356 | - info = LAPACKE_dgetrf_work(LAPACK_COL_MAJOR, sizeM, sizeN, matrix, lda, ipiv); | |
357 | - if(info != 0){ | |
358 | - stringstream ss; | |
359 | - ss << errorMessageDgetrfInfo; | |
360 | - ss << info << endl; | |
361 | - throw MolDSException(ss.str()); | |
362 | - } | |
363 | - return info; | |
364 | -} | |
365 | -} |