source: Sophya/trunk/SophyaExt/LinAlg/intflapack.cc@ 2559

Last change on this file since 2559 was 2559, checked in by cmv, 21 years ago

correction bug svddriver pour complex<r_4> et complex<r_8> (cmv 22/07/04)
File size: 30.5 KB
Line 
1#include <iostream>
2#include "intflapack.h"
3#include "tvector.h"
4#include "tmatrix.h"
5#include <typeinfo>
6
7/*************** Pour memoire (Christophe) ***************
8Les dispositions memoires (FORTRAN) pour les vecteurs et matrices LAPACK:
9
101./ --- REAL X(N):
11 if an array X of dimension (N) holds a vector x,
12 then X(i) holds "x_i" for i=1,...,N
13
142./ --- REAL A(LDA,N):
15 if a two-dimensional array A of dimension (LDA,N) holds an m-by-n matrix A,
16 then A(i,j) holds "a_ij" for i=1,...,m et j=1,...,n (LDA must be at least m).
17 Note that array arguments are usually declared in the software as assumed-size
18 arrays (last dimension *), for example: REAL A(LDA,*)
19 --- Rangement en memoire:
20 | 11 12 13 14 |
21 Ex: Real A(4,4): A = | 21 22 23 24 |
22 | 31 32 33 34 |
23 | 41 42 43 44 |
24 memoire: {11 21 31 41} {12 22 32 42} {13 23 33 43} {14 24 34 44}
25 First indice (line) "i" varies then the second (column):
26 (put all the first column, then put all the second column,
27 ..., then put all the last column)
28***********************************************************/
29
30/*!
31 \defgroup LinAlg LinAlg module
32 This module contains classes and functions for complex linear
33 algebra on arrays. This module is intended mainly to have
34 classes implementing C++ interfaces between Sophya objects
35 and external linear algebra libraries, such as LAPACK.
36*/
37
38/*!
39 \class SOPHYA::LapackServer
40 \ingroup LinAlg
41 This class implements an interface to LAPACK library driver routines.
42 The LAPACK (Linear Algebra PACKage) is a collection high performance
43 routines to solve common problems in numerical linear algebra.
44 its is available from http://www.netlib.org.
45
46 The present version of our LapackServer (Feb 2001) provides only
47 interfaces for the linear system solver and singular value
48 decomposition (SVD). Only arrays with BaseArray::FortranMemoryMapping
49 can be handled by LapackServer. LapackServer can be instanciated
50 for simple and double precision real or complex array types.
51
52 The example below shows solving a linear system A*X = B
53
54 \code
55 #include "intflapack.h"
56 // ...
57 // Use FortranMemoryMapping as default
58 BaseArray::SetDefaultMemoryMapping(BaseArray::FortranMemoryMapping);
59 // Create an fill the arrays A and B
60 int n = 20;
61 Matrix A(n, n);
62 A = RandomSequence();
63 Vector X(n),B(n);
64 X = RandomSequence();
65 B = A*X;
66 // Solve the linear system A*X = B
67 LapackServer<r_8> lps;
68 lps.LinSolve(A,B);
69 // We get the result in B, which should be equal to X ...
70 // Compute the difference B-X ;
71 Vector diff = B-X;
72 \endcode
73
74*/
75
76////////////////////////////////////////////////////////////////////////////////////
77extern "C" {
78// Le calculateur de workingspace
79 int_4 ilaenv_(int_4 *ispec,char *name,char *opts,int_4 *n1,int_4 *n2,int_4 *n3,int_4 *n4,
80 int_4 nc1,int_4 nc2);
81
82// Drivers pour resolution de systemes lineaires
83 void sgesv_(int_4* n, int_4* nrhs, r_4* a, int_4* lda,
84 int_4* ipiv, r_4* b, int_4* ldb, int_4* info);
85 void dgesv_(int_4* n, int_4* nrhs, r_8* a, int_4* lda,
86 int_4* ipiv, r_8* b, int_4* ldb, int_4* info);
87 void cgesv_(int_4* n, int_4* nrhs, complex<r_4>* a, int_4* lda,
88 int_4* ipiv, complex<r_4>* b, int_4* ldb, int_4* info);
89 void zgesv_(int_4* n, int_4* nrhs, complex<r_8>* a, int_4* lda,
90 int_4* ipiv, complex<r_8>* b, int_4* ldb, int_4* info);
91
92// Drivers pour resolution de systemes lineaires symetriques
93 void ssysv_(char* uplo, int_4* n, int_4* nrhs, r_4* a, int_4* lda,
94 int_4* ipiv, r_4* b, int_4* ldb,
95 r_4* work, int_4* lwork, int_4* info);
96 void dsysv_(char* uplo, int_4* n, int_4* nrhs, r_8* a, int_4* lda,
97 int_4* ipiv, r_8* b, int_4* ldb,
98 r_8* work, int_4* lwork, int_4* info);
99 void csysv_(char* uplo, int_4* n, int_4* nrhs, complex<r_4>* a, int_4* lda,
100 int_4* ipiv, complex<r_4>* b, int_4* ldb,
101 complex<r_4>* work, int_4* lwork, int_4* info);
102 void zsysv_(char* uplo, int_4* n, int_4* nrhs, complex<r_8>* a, int_4* lda,
103 int_4* ipiv, complex<r_8>* b, int_4* ldb,
104 complex<r_8>* work, int_4* lwork, int_4* info);
105
106// Driver pour resolution de systemes au sens de Xi2
107 void sgels_(char * trans, int_4* m, int_4* n, int_4* nrhs, r_4* a, int_4* lda,
108 r_4* b, int_4* ldb, r_4* work, int_4* lwork, int_4* info);
109 void dgels_(char * trans, int_4* m, int_4* n, int_4* nrhs, r_8* a, int_4* lda,
110 r_8* b, int_4* ldb, r_8* work, int_4* lwork, int_4* info);
111 void cgels_(char * trans, int_4* m, int_4* n, int_4* nrhs, complex<r_4>* a, int_4* lda,
112 complex<r_4>* b, int_4* ldb, complex<r_4>* work, int_4* lwork, int_4* info);
113 void zgels_(char * trans, int_4* m, int_4* n, int_4* nrhs, complex<r_8>* a, int_4* lda,
114 complex<r_8>* b, int_4* ldb, complex<r_8>* work, int_4* lwork, int_4* info);
115
116// Driver pour decomposition SVD
117 void sgesvd_(char* jobu, char* jobvt, int_4* m, int_4* n, r_4* a, int_4* lda,
118 r_4* s, r_4* u, int_4* ldu, r_4* vt, int_4* ldvt,
119 r_4* work, int_4* lwork, int_4* info);
120 void dgesvd_(char* jobu, char* jobvt, int_4* m, int_4* n, r_8* a, int_4* lda,
121 r_8* s, r_8* u, int_4* ldu, r_8* vt, int_4* ldvt,
122 r_8* work, int_4* lwork, int_4* info);
123 void cgesvd_(char* jobu, char* jobvt, int_4* m, int_4* n, complex<r_4>* a, int_4* lda,
124 r_4* s, complex<r_4>* u, int_4* ldu, complex<r_4>* vt, int_4* ldvt,
125 complex<r_4>* work, int_4* lwork, r_4* rwork, int_4* info);
126 void zgesvd_(char* jobu, char* jobvt, int_4* m, int_4* n, complex<r_8>* a, int_4* lda,
127 r_8* s, complex<r_8>* u, int_4* ldu, complex<r_8>* vt, int_4* ldvt,
128 complex<r_8>* work, int_4* lwork, r_8* rwork, int_4* info);
129
130// Driver pour eigen decomposition for symetric/hermitian matrices
131 void ssyev_(char* jobz, char* uplo, int_4* n, r_4* a, int_4* lda, r_4* w,
132 r_4* work, int_4 *lwork, int_4* info);
133 void dsyev_(char* jobz, char* uplo, int_4* n, r_8* a, int_4* lda, r_8* w,
134 r_8* work, int_4 *lwork, int_4* info);
135 void cheev_(char* jobz, char* uplo, int_4* n, complex<r_4>* a, int_4* lda, r_4* w,
136 complex<r_4>* work, int_4 *lwork, r_4* rwork, int_4* info);
137 void zheev_(char* jobz, char* uplo, int_4* n, complex<r_8>* a, int_4* lda, r_8* w,
138 complex<r_8>* work, int_4 *lwork, r_8* rwork, int_4* info);
139
140// Driver pour eigen decomposition for general squared matrices
141 void sgeev_(char* jobl, char* jobvr, int_4* n, r_4* a, int_4* lda, r_4* wr, r_4* wi,
142 r_4* vl, int_4* ldvl, r_4* vr, int_4* ldvr,
143 r_4* work, int_4 *lwork, int_4* info);
144 void dgeev_(char* jobl, char* jobvr, int_4* n, r_8* a, int_4* lda, r_8* wr, r_8* wi,
145 r_8* vl, int_4* ldvl, r_8* vr, int_4* ldvr,
146 r_8* work, int_4 *lwork, int_4* info);
147 void cgeev_(char* jobl, char* jobvr, int_4* n, complex<r_4>* a, int_4* lda, complex<r_4>* w,
148 complex<r_4>* vl, int_4* ldvl, complex<r_4>* vr, int_4* ldvr,
149 complex<r_4>* work, int_4 *lwork, r_4* rwork, int_4* info);
150 void zgeev_(char* jobl, char* jobvr, int_4* n, complex<r_8>* a, int_4* lda, complex<r_8>* w,
151 complex<r_8>* vl, int_4* ldvl, complex<r_8>* vr, int_4* ldvr,
152 complex<r_8>* work, int_4 *lwork, r_8* rwork, int_4* info);
153
154}
155
156// -------------- Classe LapackServer<T> --------------
157
158////////////////////////////////////////////////////////////////////////////////////
159template <class T>
160LapackServer<T>::LapackServer()
161{
162 SetWorkSpaceSizeFactor();
163}
164
165template <class T>
166LapackServer<T>::~LapackServer()
167{
168}
169
170// --- ATTENTION BUG POSSIBLE dans l'avenir (CMV) --- REZA A LIRE S.T.P.
171// -> Cette connerie de Fortran/C interface
172// Dans les routines fortran de lapack:
173// Appel depuis le C avec:
174// int_4 lwork = -1;
175// SUBROUTINE SSYSV( UPLO,N,NRHS,A,LDA,IPIV,B,LDB,WORK,LWORK,INFO)
176// INTEGER INFO, LDA, LDB, LWORK, N, NRHS
177// LOGICAL LQUERY
178// LQUERY = ( LWORK.EQ.-1 )
179// ELSE IF( LWORK.LT.1 .AND. .NOT.LQUERY ) THEN
180// ==> le test est bien interprete sous Linux mais pas sous OSF
181// ==> Sous OSF "LWORK.EQ.-1" est FALSE quand on passe lwork=-1 par argument
182// ==> POUR REZA: confusion entier 4 / 8 bits ??? (bizarre on l'aurait vu avant?)
183////////////////////////////////////////////////////////////////////////////////////
184template <class T>
185int_4 LapackServer<T>::ilaenv_en_C(int_4 ispec,char *name,char *opts,int_4 n1,int_4 n2,int_4 n3,int_4 n4)
186{
187 int_4 nc1 = strlen(name);
188 int_4 nc2 = strlen(opts);
189 int_4 rc=0;
190 rc = ilaenv_(&ispec,name,opts,&n1,&n2,&n3,&n4,nc1,nc2);
191 //cout<<"ilaenv_en_C("<<ispec<<","<<name<<"("<<nc1<<"),"<<opts<<"("<<nc2<<"),"
192 // <<n1<<","<<n2<<","<<n3<<","<<n4<<") = "<<rc<<endl;
193 return rc;
194}
195
196////////////////////////////////////////////////////////////////////////////////////
197//! Interface to Lapack linear system solver driver s/d/c/zgesvd().
198/*! Solve the linear system a * x = b. Input arrays
199 should have FortranMemory mapping (column packed).
200 \param a : input matrix, overwritten on output
201 \param b : input-output, input vector b, contains x on exit
202 \return : return code from lapack driver _gesv()
203 */
204template <class T>
205int LapackServer<T>::LinSolve(TArray<T>& a, TArray<T> & b)
206{
207 if ( ( a.NbDimensions() != 2 ) || ( b.NbDimensions() != 2 ) )
208 throw(SzMismatchError("LapackServer::LinSolve(a,b) a Or b NbDimensions() != 2"));
209
210 int_4 rowa = a.RowsKA();
211 int_4 cola = a.ColsKA();
212 int_4 rowb = b.RowsKA();
213 int_4 colb = b.ColsKA();
214 if ( a.Size(rowa) != a.Size(cola))
215 throw(SzMismatchError("LapackServer::LinSolve(a,b) a Not a square Array"));
216 if ( a.Size(rowa) != b.Size(rowb))
217 throw(SzMismatchError("LapackServer::LinSolve(a,b) RowSize(a <> b) "));
218
219 if (!a.IsPacked(rowa) || !b.IsPacked(rowb))
220 throw(SzMismatchError("LapackServer::LinSolve(a,b) a Or b Not Column Packed"));
221
222 int_4 n = a.Size(rowa);
223 int_4 nrhs = b.Size(colb);
224 int_4 lda = a.Step(cola);
225 int_4 ldb = b.Step(colb);
226 int_4 info;
227 int_4* ipiv = new int_4[n];
228
229 if (typeid(T) == typeid(r_4) )
230 sgesv_(&n, &nrhs, (r_4 *)a.Data(), &lda, ipiv, (r_4 *)b.Data(), &ldb, &info);
231 else if (typeid(T) == typeid(r_8) )
232 dgesv_(&n, &nrhs, (r_8 *)a.Data(), &lda, ipiv, (r_8 *)b.Data(), &ldb, &info);
233 else if (typeid(T) == typeid(complex<r_4>) )
234 cgesv_(&n, &nrhs, (complex<r_4> *)a.Data(), &lda, ipiv,
235 (complex<r_4> *)b.Data(), &ldb, &info);
236 else if (typeid(T) == typeid(complex<r_8>) )
237 zgesv_(&n, &nrhs, (complex<r_8> *)a.Data(), &lda, ipiv,
238 (complex<r_8> *)b.Data(), &ldb, &info);
239 else {
240 delete[] ipiv;
241 string tn = typeid(T).name();
242 cerr << " LapackServer::LinSolve(a,b) - Unsupported DataType T = " << tn << endl;
243 throw TypeMismatchExc("LapackServer::LinSolve(a,b) - Unsupported DataType (T)");
244 }
245 delete[] ipiv;
246 return(info);
247}
248
249////////////////////////////////////////////////////////////////////////////////////
250//! Interface to Lapack linear system solver driver s/d/c/zsysvd().
251/*! Solve the linear system a * x = b with a symetric. Input arrays
252 should have FortranMemory mapping (column packed).
253 \param a : input matrix symetric , overwritten on output
254 \param b : input-output, input vector b, contains x on exit
255 \return : return code from lapack driver _gesv()
256 */
257template <class T>
258int LapackServer<T>::LinSolveSym(TArray<T>& a, TArray<T> & b)
259// --- REMARQUES DE CMV ---
260// 1./ contrairement a ce qui est dit dans la doc, il s'agit
261// de matrices SYMETRIQUES complexes et non HERMITIENNES !!!
262// 2./ pourquoi les routines de LinSolve pour des matrices symetriques
263// sont plus de deux fois plus lentes que les LinSolve generales sur OSF
264// et sensiblement plus lentes sous Linux ???
265{
266 if ( ( a.NbDimensions() != 2 ) || ( b.NbDimensions() != 2 ) )
267 throw(SzMismatchError("LapackServer::LinSolveSym(a,b) a Or b NbDimensions() != 2"));
268 int_4 rowa = a.RowsKA();
269 int_4 cola = a.ColsKA();
270 int_4 rowb = b.RowsKA();
271 int_4 colb = b.ColsKA();
272 if ( a.Size(rowa) != a.Size(cola))
273 throw(SzMismatchError("LapackServer::LinSolveSym(a,b) a Not a square Array"));
274 if ( a.Size(rowa) != b.Size(rowb))
275 throw(SzMismatchError("LapackServer::LinSolveSym(a,b) RowSize(a <> b) "));
276
277 if (!a.IsPacked(rowa) || !b.IsPacked(rowb))
278 throw(SzMismatchError("LapackServer::LinSolveSym(a,b) a Or b Not Column Packed"));
279
280 int_4 n = a.Size(rowa);
281 int_4 nrhs = b.Size(colb);
282 int_4 lda = a.Step(cola);
283 int_4 ldb = b.Step(colb);
284 int_4 info = 0;
285 int_4* ipiv = new int_4[n];
286 int_4 lwork = -1;
287 T * work = NULL;
288
289 char uplo = 'U'; // char uplo = 'L';
290 char struplo[5]; struplo[0] = uplo; struplo[1] = '\0';
291
292 if (typeid(T) == typeid(r_4) ) {
293 lwork = ilaenv_en_C(1,"SSYTRF",struplo,n,-1,-1,-1) * n +5;
294 work = new T[lwork];
295 ssysv_(&uplo, &n, &nrhs, (r_4 *)a.Data(), &lda, ipiv, (r_4 *)b.Data(), &ldb,
296 (r_4 *)work, &lwork, &info);
297 } else if (typeid(T) == typeid(r_8) ) {
298 lwork = ilaenv_en_C(1,"DSYTRF",struplo,n,-1,-1,-1) * n +5;
299 work = new T[lwork];
300 dsysv_(&uplo, &n, &nrhs, (r_8 *)a.Data(), &lda, ipiv, (r_8 *)b.Data(), &ldb,
301 (r_8 *)work, &lwork, &info);
302 } else if (typeid(T) == typeid(complex<r_4>) ) {
303 lwork = ilaenv_en_C(1,"CSYTRF",struplo,n,-1,-1,-1) * n +5;
304 work = new T[lwork];
305 csysv_(&uplo, &n, &nrhs, (complex<r_4> *)a.Data(), &lda, ipiv,
306 (complex<r_4> *)b.Data(), &ldb,
307 (complex<r_4> *)work, &lwork, &info);
308 } else if (typeid(T) == typeid(complex<r_8>) ) {
309 lwork = ilaenv_en_C(1,"ZSYTRF",struplo,n,-1,-1,-1) * n +5;
310 work = new T[lwork];
311 zsysv_(&uplo, &n, &nrhs, (complex<r_8> *)a.Data(), &lda, ipiv,
312 (complex<r_8> *)b.Data(), &ldb,
313 (complex<r_8> *)work, &lwork, &info);
314 } else {
315 if(work) delete[] work;
316 delete[] ipiv;
317 string tn = typeid(T).name();
318 cerr << " LapackServer::LinSolveSym(a,b) - Unsupported DataType T = " << tn << endl;
319 throw TypeMismatchExc("LapackServer::LinSolveSym(a,b) - Unsupported DataType (T)");
320 }
321 if(work) delete[] work;
322 delete[] ipiv;
323 return(info);
324}
325
326////////////////////////////////////////////////////////////////////////////////////
327//! Interface to Lapack least squares solver driver s/d/c/zgels().
328/*! Solves the linear least squares problem defined by an m-by-n matrix
329 \b a and an m element vector \b b .
330 A solution \b x to the overdetermined system of linear equations
331 b = a * x is computed, minimizing the norm of b-a*x.
332 Underdetermined systems (m<n) are not yet handled.
333 Inout arrays should have FortranMemory mapping (column packed).
334 \param a : input matrix, overwritten on output
335 \param b : input-output, input vector b, contains x on exit.
336 \return : return code from lapack driver _gels()
337 \warning : b is not resized.
338 */
339/*
340 $CHECK$ - A faire - cas m<n
341 If the linear system is underdetermined, the minimum norm
342 solution is computed.
343*/
344
345template <class T>
346int LapackServer<T>::LeastSquareSolve(TArray<T>& a, TArray<T> & b)
347{
348 if ( ( a.NbDimensions() != 2 ) || ( b.NbDimensions() != 2 ) )
349 throw(SzMismatchError("LapackServer::LinSolve(a,b) a Or b NbDimensions() != 2"));
350
351 int_4 rowa = a.RowsKA();
352 int_4 cola = a.ColsKA();
353 int_4 rowb = b.RowsKA();
354 int_4 colb = b.ColsKA();
355
356
357 if ( a.Size(rowa) != b.Size(rowb))
358 throw(SzMismatchError("LapackServer::LeastSquareSolve(a,b) RowSize(a <> b) "));
359
360 if (!a.IsPacked(rowa) || !b.IsPacked(rowb))
361 throw(SzMismatchError("LapackServer::LeastSquareSolve(a,b) a Or b Not Column Packed"));
362
363 if ( a.Size(rowa) < a.Size(cola)) { // $CHECK$ - m<n a changer
364 cout << " LapackServer<T>::LeastSquareSolve() - m<n - Not yet implemented for "
365 << " underdetermined systems ! " << endl;
366 throw(SzMismatchError("LapackServer::LeastSquareSolve(a,b) NRows<NCols - "));
367 }
368 int_4 m = a.Size(rowa);
369 int_4 n = a.Size(cola);
370 int_4 nrhs = b.Size(colb);
371
372 int_4 lda = a.Step(cola);
373 int_4 ldb = b.Step(colb);
374 int_4 info;
375
376 int_4 minmn = (m < n) ? m : n;
377 int_4 maxmn = (m > n) ? m : n;
378 int_4 maxmnrhs = (nrhs > maxmn) ? nrhs : maxmn;
379 if (maxmnrhs < 1) maxmnrhs = 1;
380
381 int_4 lwork = minmn+maxmnrhs*5;
382 T * work = new T[lwork];
383
384 char trans = 'N';
385
386 if (typeid(T) == typeid(r_4) )
387 sgels_(&trans, &m, &n, &nrhs, (r_4 *)a.Data(), &lda,
388 (r_4 *)b.Data(), &ldb, (r_4 *)work, &lwork, &info);
389 else if (typeid(T) == typeid(r_8) )
390 dgels_(&trans, &m, &n, &nrhs, (r_8 *)a.Data(), &lda,
391 (r_8 *)b.Data(), &ldb, (r_8 *)work, &lwork, &info);
392 else if (typeid(T) == typeid(complex<r_4>) )
393 cgels_(&trans, &m, &n, &nrhs, (complex<r_4> *)a.Data(), &lda,
394 (complex<r_4> *)b.Data(), &ldb, (complex<r_4> *)work, &lwork, &info);
395 else if (typeid(T) == typeid(complex<r_8>) )
396 zgels_(&trans, &m, &n, &nrhs, (complex<r_8> *)a.Data(), &lda,
397 (complex<r_8> *)b.Data(), &ldb, (complex<r_8> *)work, &lwork, &info);
398 else {
399 delete[] work;
400 string tn = typeid(T).name();
401 cerr << " LapackServer::LeastSquareSolve(a,b) - Unsupported DataType T = " << tn << endl;
402 throw TypeMismatchExc("LapackServer::LeastSquareSolve(a,b) - Unsupported DataType (T)");
403 }
404 delete[] work;
405 return(info);
406}
407
408
409////////////////////////////////////////////////////////////////////////////////////
410//! Interface to Lapack SVD driver s/d/c/zgesv().
411/*! Computes the vector of singular values of \b a. Input arrays
412 should have FortranMemoryMapping (column packed).
413 \param a : input m-by-n matrix
414 \param s : Vector of min(m,n) singular values (descending order)
415 \return : return code from lapack driver _gesvd()
416 */
417
418template <class T>
419int LapackServer<T>::SVD(TArray<T>& a, TArray<T> & s)
420{
421 return (SVDDriver(a, s, NULL, NULL) );
422}
423
424//! Interface to Lapack SVD driver s/d/c/zgesv().
425/*! Computes the vector of singular values of \b a, as well as
426 right and left singular vectors of \b a.
427 \f[
428 A = U \Sigma V^T , ( A = U \Sigma V^H \ complex)
429 \f]
430 \f[
431 A v_i = \sigma_i u_i \ and A^T u_i = \sigma_i v_i \ (A^H \ complex)
432 \f]
433 U and V are orthogonal (unitary) matrices.
434 \param a : input m-by-n matrix (in FotranMemoryMapping)
435 \param s : Vector of min(m,n) singular values (descending order)
436 \param u : Matrix of left singular vectors
437 \param vt : Transpose of right singular vectors.
438 \return : return code from lapack driver _gesvd()
439 */
440template <class T>
441int LapackServer<T>::SVD(TArray<T>& a, TArray<T> & s, TArray<T> & u, TArray<T> & vt)
442{
443 return (SVDDriver(a, s, &u, &vt) );
444}
445
446
447//! Interface to Lapack SVD driver s/d/c/zgesv().
448template <class T>
449int LapackServer<T>::SVDDriver(TArray<T>& a, TArray<T> & s, TArray<T>* up, TArray<T>* vtp)
450{
451 if ( ( a.NbDimensions() != 2 ) )
452 throw(SzMismatchError("LapackServer::SVD(a, ...) a.NbDimensions() != 2"));
453
454 int_4 rowa = a.RowsKA();
455 int_4 cola = a.ColsKA();
456
457 if ( !a.IsPacked(rowa) )
458 throw(SzMismatchError("LapackServer::SVD(a, ...) a Not Column Packed "));
459
460 int_4 m = a.Size(rowa);
461 int_4 n = a.Size(cola);
462 int_4 maxmn = (m > n) ? m : n;
463 int_4 minmn = (m < n) ? m : n;
464
465 char jobu, jobvt;
466 jobu = 'N';
467 jobvt = 'N';
468
469 sa_size_t sz[2];
470 if ( up != NULL) {
471 if ( dynamic_cast< TVector<T> * > (vtp) )
472 throw( TypeMismatchExc("LapackServer::SVD() Wrong type (=TVector<T>) for u !") );
473 up->SetMemoryMapping(BaseArray::FortranMemoryMapping);
474 sz[0] = sz[1] = m;
475 up->ReSize(2, sz );
476 jobu = 'A';
477 }
478 else {
479 up = new TMatrix<T>(1,1);
480 jobu = 'N';
481 }
482 if ( vtp != NULL) {
483 if ( dynamic_cast< TVector<T> * > (vtp) )
484 throw( TypeMismatchExc("LapackServer::SVD() Wrong type (=TVector<T>) for vt !") );
485 vtp->SetMemoryMapping(BaseArray::FortranMemoryMapping);
486 sz[0] = sz[1] = n;
487 vtp->ReSize(2, sz );
488 jobvt = 'A';
489 }
490 else {
491 vtp = new TMatrix<T>(1,1);
492 jobvt = 'N';
493 }
494
495 TVector<T> *vs = dynamic_cast< TVector<T> * > (&s);
496 if (vs) vs->ReSize(minmn);
497 else {
498 TMatrix<T> *ms = dynamic_cast< TMatrix<T> * > (&s);
499 if (ms) ms->ReSize(minmn,1);
500 else {
501 sz[0] = minmn; sz[1] = 1;
502 s.ReSize(1, sz);
503 }
504 }
505
506 int_4 lda = a.Step(a.ColsKA());
507 int_4 ldu = up->Step(up->ColsKA());
508 int_4 ldvt = vtp->Step(vtp->ColsKA());
509
510 int_4 lwork = maxmn*5*wspace_size_factor;
511 T * work = new T[lwork];
512 int_4 info;
513
514 if (typeid(T) == typeid(r_4) ) {
515 sgesvd_(&jobu, &jobvt, &m, &n, (r_4 *)a.Data(), &lda,
516 (r_4 *)s.Data(), (r_4 *) up->Data(), &ldu, (r_4 *)vtp->Data(), &ldvt,
517 (r_4 *)work, &lwork, &info);
518 } else if (typeid(T) == typeid(r_8) ) {
519 dgesvd_(&jobu, &jobvt, &m, &n, (r_8 *)a.Data(), &lda,
520 (r_8 *)s.Data(), (r_8 *) up->Data(), &ldu, (r_8 *)vtp->Data(), &ldvt,
521 (r_8 *)work, &lwork, &info);
522 } else if (typeid(T) == typeid(complex<r_4>) ) {
523 r_4 * rwork = new r_4[5*minmn +5];
524 r_4 * sloc = new r_4[minmn];
525 cgesvd_(&jobu, &jobvt, &m, &n, (complex<r_4> *)a.Data(), &lda,
526 (r_4 *)sloc, (complex<r_4> *) up->Data(), &ldu,
527 (complex<r_4> *)vtp->Data(), &ldvt,
528 (complex<r_4> *)work, &lwork, (r_4 *)rwork, &info);
529 for(int_4 i=0;i<minmn;i++) s[i] = sloc[i];
530 delete [] rwork; delete [] sloc;
531 } else if (typeid(T) == typeid(complex<r_8>) ) {
532 r_8 * rwork = new r_8[5*minmn +5];
533 r_8 * sloc = new r_8[minmn];
534 zgesvd_(&jobu, &jobvt, &m, &n, (complex<r_8> *)a.Data(), &lda,
535 (r_8 *)sloc, (complex<r_8> *) up->Data(), &ldu,
536 (complex<r_8> *)vtp->Data(), &ldvt,
537 (complex<r_8> *)work, &lwork, (r_8 *)rwork, &info);
538 for(int_4 i=0;i<minmn;i++) s[i] = sloc[i];
539 delete [] rwork; delete [] sloc;
540 } else {
541 if (jobu == 'N') delete up;
542 if (jobvt == 'N') delete vtp;
543 string tn = typeid(T).name();
544 cerr << " LapackServer::SVDDriver(...) - Unsupported DataType T = " << tn << endl;
545 throw TypeMismatchExc("LapackServer::LinSolve(a,b) - Unsupported DataType (T)");
546 }
547
548 if (jobu == 'N') delete up;
549 if (jobvt == 'N') delete vtp;
550 return(info);
551}
552
553
554////////////////////////////////////////////////////////////////////////////////////
555/*! Computes the eigen values and eigen vectors of a symetric (or hermitian) matrix \b a.
556 Input arrays should have FortranMemoryMapping (column packed).
557 \param a : input symetric (or hermitian) n-by-n matrix
558 \param b : Vector of eigenvalues (descending order)
559 \param eigenvector : if true compute eigenvectors, if not only eigenvalues
560 \param a : on return array of eigenvectors (same order than eval, one vector = one column)
561 \return : return code from lapack driver _gesvd()
562 */
563
564template <class T>
565int LapackServer<T>::LapackEigenSym(TArray<T>& a, TVector<r_8>& b, bool eigenvector)
566{
567 if ( a.NbDimensions() != 2 )
568 throw(SzMismatchError("LapackServer::LapackEigenSym(a,b) a NbDimensions() != 2"));
569 int_4 rowa = a.RowsKA();
570 int_4 cola = a.ColsKA();
571 if ( a.Size(rowa) != a.Size(cola))
572 throw(SzMismatchError("LapackServer::LapackEigenSym(a,b) a Not a square Array"));
573 if (!a.IsPacked(rowa))
574 throw(SzMismatchError("LapackServer::LapackEigenSym(a,b) a Not Column Packed"));
575
576 char uplo='U'; char struplo[5]; struplo[0]=uplo; struplo[1]='\0';
577 char jobz='N'; if(eigenvector) jobz='V';
578 char strjobz[5]; strjobz[0]=jobz; strjobz[1]='\0';
579
580 int_4 n = a.Size(rowa);
581 int_4 lda = a.Step(cola);
582 int_4 info = 0;
583
584 b.ReSize(n); b = 0.;
585
586 if (typeid(T) == typeid(r_4) ) {
587 int_4 lwork = 3*n-1 +5; r_4* work = new r_4[lwork];
588 r_4* w = new r_4[n];
589 ssyev_(strjobz,struplo,&n,(r_4 *)a.Data(),&lda,(r_4 *)w,(r_4 *)work,&lwork,&info);
590 if(info==0) for(int i=0;i<n;i++) b(i) = w[i];
591 delete [] work; delete [] w;
592 } else if (typeid(T) == typeid(r_8) ) {
593 int_4 lwork = 3*n-1 +5; r_8* work = new r_8[lwork];
594 r_8* w = new r_8[n];
595 dsyev_(strjobz,struplo,&n,(r_8 *)a.Data(),&lda,(r_8 *)w,(r_8 *)work,&lwork,&info);
596 if(info==0) for(int i=0;i<n;i++) b(i) = w[i];
597 delete [] work; delete [] w;
598 } else if (typeid(T) == typeid(complex<r_4>) ) {
599 int_4 lwork = 2*n-1 +5; complex<r_4>* work = new complex<r_4>[lwork];
600 r_4* rwork = new r_4[3*n-2 +5]; r_4* w = new r_4[n];
601 cheev_(strjobz,struplo,&n,(complex<r_4> *)a.Data(),&lda,(r_4 *)w
602 ,(complex<r_4> *)work,&lwork,(r_4 *)rwork,&info);
603 if(info==0) for(int i=0;i<n;i++) b(i) = w[i];
604 delete [] work; delete [] rwork; delete [] w;
605 } else if (typeid(T) == typeid(complex<r_8>) ) {
606 int_4 lwork = 2*n-1 +5; complex<r_8>* work = new complex<r_8>[lwork];
607 r_8* rwork = new r_8[3*n-2 +5]; r_8* w = new r_8[n];
608 zheev_(strjobz,struplo,&n,(complex<r_8> *)a.Data(),&lda,(r_8 *)w
609 ,(complex<r_8> *)work,&lwork,(r_8 *)rwork,&info);
610 if(info==0) for(int i=0;i<n;i++) b(i) = w[i];
611 delete [] work; delete [] rwork; delete [] w;
612 } else {
613 string tn = typeid(T).name();
614 cerr << " LapackServer::LapackEigenSym(a,b) - Unsupported DataType T = " << tn << endl;
615 throw TypeMismatchExc("LapackServer::LapackEigenSym(a,b) - Unsupported DataType (T)");
616 }
617
618 return(info);
619}
620
621////////////////////////////////////////////////////////////////////////////////////
622/*! Computes the eigen values and eigen vectors of a general squared matrix \b a.
623 Input arrays should have FortranMemoryMapping (column packed).
624 \param a : input general n-by-n matrix
625 \param eval : Vector of eigenvalues (complex double precision)
626 \param evec : Matrix of eigenvector (same order than eval, one vector = one column)
627 \param eigenvector : if true compute (right) eigenvectors, if not only eigenvalues
628 \param a : on return array of eigenvectors
629 \return : return code from lapack driver _gesvd()
630 \verbatim
631 eval : contains the computed eigenvalues.
632 --- For real matrices "a" :
633 Complex conjugate pairs of eigenvalues appear consecutively
634 with the eigenvalue having the positive imaginary part first.
635 evec : the right eigenvectors v(j) are stored one after another
636 in the columns of evec, in the same order as their eigenvalues.
637 --- For real matrices "a" :
638 If the j-th eigenvalue is real, then v(j) = evec(:,j),
639 the j-th column of evec.
640 If the j-th and (j+1)-st eigenvalues form a complex
641 conjugate pair, then v(j) = evec(:,j) + i*evec(:,j+1) and
642 v(j+1) = evec(:,j) - i*evec(:,j+1).
643 \endverbatim
644*/
645
646template <class T>
647int LapackServer<T>::LapackEigen(TArray<T>& a, TVector< complex<r_8> >& eval, TMatrix<T>& evec, bool eigenvector)
648{
649 if ( a.NbDimensions() != 2 )
650 throw(SzMismatchError("LapackServer::LapackEigen(a,b) a NbDimensions() != 2"));
651 int_4 rowa = a.RowsKA();
652 int_4 cola = a.ColsKA();
653 if ( a.Size(rowa) != a.Size(cola))
654 throw(SzMismatchError("LapackServer::LapackEigen(a,b) a Not a square Array"));
655 if (!a.IsPacked(rowa))
656 throw(SzMismatchError("LapackServer::LapackEigen(a,b) a Not Column Packed"));
657
658 char jobvl = 'N'; char strjobvl[5]; strjobvl[0] = jobvl; strjobvl[1] = '\0';
659 char jobvr = 'N'; if(eigenvector) jobvr='V';
660 char strjobvr[5]; strjobvr[0] = jobvr; strjobvr[1] = '\0';
661
662 int_4 n = a.Size(rowa);
663 int_4 lda = a.Step(cola);
664 int_4 info = 0;
665
666 eval.ReSize(n); eval = complex<r_8>(0.,0.);
667 if(eigenvector) {evec.ReSize(n,n); evec = (T) 0.;}
668 int_4 ldvr = n, ldvl = 1;
669
670 if (typeid(T) == typeid(r_4) ) {
671 int_4 lwork = 4*n +5; r_4* work = new r_4[lwork];
672 r_4* wr = new r_4[n]; r_4* wi = new r_4[n]; r_4* vl = NULL;
673 sgeev_(strjobvl,strjobvr,&n,(r_4 *)a.Data(),&lda,(r_4 *)wr,(r_4 *)wi,
674 (r_4 *)vl,&ldvl,(r_4 *)evec.Data(),&ldvr,
675 (r_4 *)work,&lwork,&info);
676 if(info==0) for(int i=0;i<n;i++) eval(i) = complex<r_8>(wr[i],wi[i]);
677 delete [] work; delete [] wr; delete [] wi;
678 } else if (typeid(T) == typeid(r_8) ) {
679 int_4 lwork = 4*n +5; r_8* work = new r_8[lwork];
680 r_8* wr = new r_8[n]; r_8* wi = new r_8[n]; r_8* vl = NULL;
681 dgeev_(strjobvl,strjobvr,&n,(r_8 *)a.Data(),&lda,(r_8 *)wr,(r_8 *)wi,
682 (r_8 *)vl,&ldvl,(r_8 *)evec.Data(),&ldvr,
683 (r_8 *)work,&lwork,&info);
684 if(info==0) for(int i=0;i<n;i++) eval(i) = complex<r_8>(wr[i],wi[i]);
685 delete [] work; delete [] wr; delete [] wi;
686 } else if (typeid(T) == typeid(complex<r_4>) ) {
687 int_4 lwork = 2*n +5; complex<r_4>* work = new complex<r_4>[lwork];
688 r_4* rwork = new r_4[2*n+5]; r_4* vl = NULL; TVector< complex<r_4> > w(n);
689 cgeev_(strjobvl,strjobvr,&n,(complex<r_4> *)a.Data(),&lda,(complex<r_4> *)w.Data(),
690 (complex<r_4> *)vl,&ldvl,(complex<r_4> *)evec.Data(),&ldvr,
691 (complex<r_4> *)work,&lwork,(r_4 *)rwork,&info);
692 if(info==0) for(int i=0;i<n;i++) eval(i) = w(i);
693 delete [] work; delete [] rwork;
694 } else if (typeid(T) == typeid(complex<r_8>) ) {
695 int_4 lwork = 2*n +5; complex<r_8>* work = new complex<r_8>[lwork];
696 r_8* rwork = new r_8[2*n+5]; r_8* vl = NULL;
697 zgeev_(strjobvl,strjobvr,&n,(complex<r_8> *)a.Data(),&lda,(complex<r_8> *)eval.Data(),
698 (complex<r_8> *)vl,&ldvl,(complex<r_8> *)evec.Data(),&ldvr,
699 (complex<r_8> *)work,&lwork,(r_8 *)rwork,&info);
700 delete [] work; delete [] rwork;
701 } else {
702 string tn = typeid(T).name();
703 cerr << " LapackServer::LapackEigen(a,b) - Unsupported DataType T = " << tn << endl;
704 throw TypeMismatchExc("LapackServer::LapackEigen(a,b) - Unsupported DataType (T)");
705 }
706
707 return(info);
708}
709
710
711
712
713
714////////////////////////////////////////////////////////////////////////////////////
715void rztest_lapack(TArray<r_4>& aa, TArray<r_4>& bb)
716{
717 if ( aa.NbDimensions() != 2 ) throw(SzMismatchError("rztest_lapack(TMatrix<r_4> A Not a Matrix"));
718 if ( aa.SizeX() != aa.SizeY()) throw(SzMismatchError("rztest_lapack(TMatrix<r_4> A Not a square Matrix"));
719 if ( bb.NbDimensions() != 2 ) throw(SzMismatchError("rztest_lapack(TMatrix<r_4> A Not a Matrix"));
720 if ( bb.SizeX() != aa.SizeX() ) throw(SzMismatchError("rztest_lapack(TMatrix<r_4> A <> B "));
721 if ( !bb.IsPacked() || !bb.IsPacked() )
722 throw(SzMismatchError("rztest_lapack(TMatrix<r_4> Not packed A or B "));
723
724 int_4 n = aa.SizeX();
725 int_4 nrhs = bb.SizeY();
726 int_4 lda = n;
727 int_4 ldb = bb.SizeX();
728 int_4 info;
729 int_4* ipiv = new int_4[n];
730 sgesv_(&n, &nrhs, aa.Data(), &lda, ipiv, bb.Data(), &ldb, &info);
731 delete[] ipiv;
732 cout << "rztest_lapack/Info= " << info << endl;
733 cout << aa << "\n" << bb << endl;
734 return;
735}
736
737///////////////////////////////////////////////////////////////
738#ifdef __CXX_PRAGMA_TEMPLATES__
739#pragma define_template LapackServer<r_4>
740#pragma define_template LapackServer<r_8>
741#pragma define_template LapackServer< complex<r_4> >
742#pragma define_template LapackServer< complex<r_8> >
743#endif
744
745#if defined(ANSI_TEMPLATES) || defined(GNU_TEMPLATES)
746template class LapackServer<r_4>;
747template class LapackServer<r_8>;
748template class LapackServer< complex<r_4> >;
749template class LapackServer< complex<r_8> >;
750#endif
751
752#if defined(OS_LINUX)
753// Pour le link avec f2c sous Linux
754extern "C" {
755 void MAIN__();
756}
757
758void MAIN__()
759{
760 cerr << "MAIN__() function for linking with libf2c.a " << endl;
761 cerr << " This function should never be called !!! " << endl;
762 throw PError("MAIN__() should not be called - see intflapack.cc");
763}
764#endif
Note: See TracBrowser for help on using the repository browser.