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

Last change on this file since 2864 was 2646, checked in by ansari, 21 years ago

Ajout methode calcul de la matrice inverse par Lapack en utilisant la resolution de systeme + MAJ doc - Reza 7 Fev 2005
File size: 45.5 KB
Line 
1#include <iostream>
2#include <math.h>
3#include "sopnamsp.h"
4#include "intflapack.h"
5#include "tvector.h"
6#include "tmatrix.h"
7#include <typeinfo>
8
9#define GARDMEM 5
10
11/*************** Pour memoire (Christophe) ***************
12Les dispositions memoires (FORTRAN) pour les vecteurs et matrices LAPACK:
13
141./ --- REAL X(N):
15 if an array X of dimension (N) holds a vector x,
16 then X(i) holds "x_i" for i=1,...,N
17
182./ --- REAL A(LDA,N):
19 if a two-dimensional array A of dimension (LDA,N) holds an m-by-n matrix A,
20 then A(i,j) holds "a_ij" for i=1,...,m et j=1,...,n (LDA must be at least m).
21 Note that array arguments are usually declared in the software as assumed-size
22 arrays (last dimension *), for example: REAL A(LDA,*)
23 --- Rangement en memoire:
24 | 11 12 13 14 |
25 Ex: Real A(4,4): A = | 21 22 23 24 |
26 | 31 32 33 34 |
27 | 41 42 43 44 |
28 memoire: {11 21 31 41} {12 22 32 42} {13 23 33 43} {14 24 34 44}
29 First indice (line) "i" varies then the second (column):
30 (put all the first column, then put all the second column,
31 ..., then put all the last column)
32***********************************************************/
33
34/*!
35 \defgroup LinAlg LinAlg module
36 This module contains classes and functions for complex linear
37 algebra on arrays. This module is intended mainly to have
38 classes implementing C++ interfaces between Sophya objects
39 and external linear algebra libraries, such as LAPACK.
40*/
41
42/*!
43 \class SOPHYA::LapackServer
44 \ingroup LinAlg
45 This class implements an interface to LAPACK library driver routines.
46 The LAPACK (Linear Algebra PACKage) is a collection high performance
47 routines to solve common problems in numerical linear algebra.
48 its is available from http://www.netlib.org.
49
50 The present version of LapackServer (Feb 2005) provides
51 interfaces for the linear system solver, singular value
52 decomposition (SVD), Least square solver and
53 eigen value / eigen vector decomposition.
54 Only arrays with BaseArray::FortranMemoryMapping
55 can be handled by LapackServer. LapackServer can be instanciated
56 for simple and double precision real or complex array types.
57 \warning The input array is overwritten in most cases.
58 The example below shows solving a linear system A*X = B
59
60 \code
61 #include "intflapack.h"
62 // ...
63 // Use FortranMemoryMapping as default
64 BaseArray::SetDefaultMemoryMapping(BaseArray::FortranMemoryMapping);
65 // Create an fill the arrays A and B
66 int n = 20;
67 Matrix A(n, n);
68 A = RandomSequence();
69 Vector X(n),B(n);
70 X = RandomSequence();
71 B = A*X;
72 // Solve the linear system A*X = B
73 LapackServer<r_8> lps;
74 lps.LinSolve(A,B);
75 // We get the result in B, which should be equal to X ...
76 // Compute the difference B-X ;
77 Vector diff = B-X;
78 \endcode
79
80*/
81
82////////////////////////////////////////////////////////////////////////////////////
83extern "C" {
84// Le calculateur de workingspace
85 int_4 ilaenv_(int_4 *ispec,char *name,char *opts,int_4 *n1,int_4 *n2,int_4 *n3,int_4 *n4,
86 int_4 nc1,int_4 nc2);
87
88// Drivers pour resolution de systemes lineaires
89 void sgesv_(int_4* n, int_4* nrhs, r_4* a, int_4* lda,
90 int_4* ipiv, r_4* b, int_4* ldb, int_4* info);
91 void dgesv_(int_4* n, int_4* nrhs, r_8* a, int_4* lda,
92 int_4* ipiv, r_8* b, int_4* ldb, int_4* info);
93 void cgesv_(int_4* n, int_4* nrhs, complex<r_4>* a, int_4* lda,
94 int_4* ipiv, complex<r_4>* b, int_4* ldb, int_4* info);
95 void zgesv_(int_4* n, int_4* nrhs, complex<r_8>* a, int_4* lda,
96 int_4* ipiv, complex<r_8>* b, int_4* ldb, int_4* info);
97
98// Drivers pour resolution de systemes lineaires symetriques
99 void ssysv_(char* uplo, int_4* n, int_4* nrhs, r_4* a, int_4* lda,
100 int_4* ipiv, r_4* b, int_4* ldb,
101 r_4* work, int_4* lwork, int_4* info);
102 void dsysv_(char* uplo, int_4* n, int_4* nrhs, r_8* a, int_4* lda,
103 int_4* ipiv, r_8* b, int_4* ldb,
104 r_8* work, int_4* lwork, int_4* info);
105 void csysv_(char* uplo, int_4* n, int_4* nrhs, complex<r_4>* a, int_4* lda,
106 int_4* ipiv, complex<r_4>* b, int_4* ldb,
107 complex<r_4>* work, int_4* lwork, int_4* info);
108 void zsysv_(char* uplo, int_4* n, int_4* nrhs, complex<r_8>* a, int_4* lda,
109 int_4* ipiv, complex<r_8>* b, int_4* ldb,
110 complex<r_8>* work, int_4* lwork, int_4* info);
111
112// Driver pour resolution de systemes au sens de Xi2
113 void sgels_(char * trans, int_4* m, int_4* n, int_4* nrhs, r_4* a, int_4* lda,
114 r_4* b, int_4* ldb, r_4* work, int_4* lwork, int_4* info);
115 void dgels_(char * trans, int_4* m, int_4* n, int_4* nrhs, r_8* a, int_4* lda,
116 r_8* b, int_4* ldb, r_8* work, int_4* lwork, int_4* info);
117 void cgels_(char * trans, int_4* m, int_4* n, int_4* nrhs, complex<r_4>* a, int_4* lda,
118 complex<r_4>* b, int_4* ldb, complex<r_4>* work, int_4* lwork, int_4* info);
119 void zgels_(char * trans, int_4* m, int_4* n, int_4* nrhs, complex<r_8>* a, int_4* lda,
120 complex<r_8>* b, int_4* ldb, complex<r_8>* work, int_4* lwork, int_4* info);
121
122// Driver pour resolution de systemes au sens de Xi2 par SVD Divide & Conquer
123 void sgelsd_(int_4* m,int_4* n,int_4* nrhs,r_4* a,int_4* lda,
124 r_4* b,int_4* ldb,r_4* s,r_4* rcond,int_4* rank,
125 r_4* work,int_4* lwork,int_4* iwork,int_4* info);
126 void dgelsd_(int_4* m,int_4* n,int_4* nrhs,r_8* a,int_4* lda,
127 r_8* b,int_4* ldb,r_8* s,r_8* rcond,int_4* rank,
128 r_8* work,int_4* lwork,int_4* iwork,int_4* info);
129 void cgelsd_(int_4* m,int_4* n,int_4* nrhs,complex<r_4>* a,int_4* lda,
130 complex<r_4>* b,int_4* ldb,r_4* s,r_4* rcond,int_4* rank,
131 complex<r_4>* work,int_4* lwork,r_4* rwork,int_4* iwork,int_4* info);
132 void zgelsd_(int_4* m,int_4* n,int_4* nrhs,complex<r_8>* a,int_4* lda,
133 complex<r_8>* b,int_4* ldb,r_8* s,r_8* rcond,int_4* rank,
134 complex<r_8>* work,int_4* lwork,r_8* rwork,int_4* iwork,int_4* info);
135
136// Driver pour decomposition SVD
137 void sgesvd_(char* jobu, char* jobvt, int_4* m, int_4* n, r_4* a, int_4* lda,
138 r_4* s, r_4* u, int_4* ldu, r_4* vt, int_4* ldvt,
139 r_4* work, int_4* lwork, int_4* info);
140 void dgesvd_(char* jobu, char* jobvt, int_4* m, int_4* n, r_8* a, int_4* lda,
141 r_8* s, r_8* u, int_4* ldu, r_8* vt, int_4* ldvt,
142 r_8* work, int_4* lwork, int_4* info);
143 void cgesvd_(char* jobu, char* jobvt, int_4* m, int_4* n, complex<r_4>* a, int_4* lda,
144 r_4* s, complex<r_4>* u, int_4* ldu, complex<r_4>* vt, int_4* ldvt,
145 complex<r_4>* work, int_4* lwork, r_4* rwork, int_4* info);
146 void zgesvd_(char* jobu, char* jobvt, int_4* m, int_4* n, complex<r_8>* a, int_4* lda,
147 r_8* s, complex<r_8>* u, int_4* ldu, complex<r_8>* vt, int_4* ldvt,
148 complex<r_8>* work, int_4* lwork, r_8* rwork, int_4* info);
149
150// Driver pour decomposition SVD Divide and Conquer
151 void sgesdd_(char* jobz, int_4* m, int_4* n, r_4* a, int_4* lda,
152 r_4* s, r_4* u, int_4* ldu, r_4* vt, int_4* ldvt,
153 r_4* work, int_4* lwork, int_4* iwork, int_4* info);
154 void dgesdd_(char* jobz, int_4* m, int_4* n, r_8* a, int_4* lda,
155 r_8* s, r_8* u, int_4* ldu, r_8* vt, int_4* ldvt,
156 r_8* work, int_4* lwork, int_4* iwork, int_4* info);
157 void cgesdd_(char* jobz, int_4* m, int_4* n, complex<r_4>* a, int_4* lda,
158 r_4* s, complex<r_4>* u, int_4* ldu, complex<r_4>* vt, int_4* ldvt,
159 complex<r_4>* work, int_4* lwork, r_4* rwork, int_4* iwork, int_4* info);
160 void zgesdd_(char* jobz, int_4* m, int_4* n, complex<r_8>* a, int_4* lda,
161 r_8* s, complex<r_8>* u, int_4* ldu, complex<r_8>* vt, int_4* ldvt,
162 complex<r_8>* work, int_4* lwork, r_8* rwork, int_4* iwork, int_4* info);
163
164// Driver pour eigen decomposition for symetric/hermitian matrices
165 void ssyev_(char* jobz, char* uplo, int_4* n, r_4* a, int_4* lda, r_4* w,
166 r_4* work, int_4 *lwork, int_4* info);
167 void dsyev_(char* jobz, char* uplo, int_4* n, r_8* a, int_4* lda, r_8* w,
168 r_8* work, int_4 *lwork, int_4* info);
169 void cheev_(char* jobz, char* uplo, int_4* n, complex<r_4>* a, int_4* lda, r_4* w,
170 complex<r_4>* work, int_4 *lwork, r_4* rwork, int_4* info);
171 void zheev_(char* jobz, char* uplo, int_4* n, complex<r_8>* a, int_4* lda, r_8* w,
172 complex<r_8>* work, int_4 *lwork, r_8* rwork, int_4* info);
173
174// Driver pour eigen decomposition for general squared matrices
175 void sgeev_(char* jobl, char* jobvr, int_4* n, r_4* a, int_4* lda, r_4* wr, r_4* wi,
176 r_4* vl, int_4* ldvl, r_4* vr, int_4* ldvr,
177 r_4* work, int_4 *lwork, int_4* info);
178 void dgeev_(char* jobl, char* jobvr, int_4* n, r_8* a, int_4* lda, r_8* wr, r_8* wi,
179 r_8* vl, int_4* ldvl, r_8* vr, int_4* ldvr,
180 r_8* work, int_4 *lwork, int_4* info);
181 void cgeev_(char* jobl, char* jobvr, int_4* n, complex<r_4>* a, int_4* lda, complex<r_4>* w,
182 complex<r_4>* vl, int_4* ldvl, complex<r_4>* vr, int_4* ldvr,
183 complex<r_4>* work, int_4 *lwork, r_4* rwork, int_4* info);
184 void zgeev_(char* jobl, char* jobvr, int_4* n, complex<r_8>* a, int_4* lda, complex<r_8>* w,
185 complex<r_8>* vl, int_4* ldvl, complex<r_8>* vr, int_4* ldvr,
186 complex<r_8>* work, int_4 *lwork, r_8* rwork, int_4* info);
187
188}
189
190// -------------- Classe LapackServer<T> --------------
191
192////////////////////////////////////////////////////////////////////////////////////
193template <class T>
194LapackServer<T>::LapackServer()
195{
196 SetWorkSpaceSizeFactor();
197}
198
199template <class T>
200LapackServer<T>::~LapackServer()
201{
202}
203
204////////////////////////////////////////////////////////////////////////////////////
205template <class T>
206int_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)
207{
208 int_4 nc1 = strlen(name), nc2 = strlen(opts), rc=0;
209 rc = ilaenv_(&ispec,name,opts,&n1,&n2,&n3,&n4,nc1,nc2);
210 //cout<<"ilaenv_en_C("<<ispec<<","<<name<<"("<<nc1<<"),"<<opts<<"("<<nc2<<"),"
211 // <<n1<<","<<n2<<","<<n3<<","<<n4<<") = "<<rc<<endl;
212 return rc;
213}
214
215template <class T>
216int_4 LapackServer<T>::type2i4(void *val,int nbytes)
217// Retourne un entier contenant la valeur contenue dans val
218// - nbytes = nombre de bytes dans le contenu de val
219// ex: r_4 x = 3.4; type2i4(&x,4) -> 3
220// ex: r_8 x = 3.4; type2i4(&x,8) -> 3
221// ex: complex<r_4> x(3.4,7.8); type2i4(&x,4) -> 3
222// ex: complex<r_8> x(3.4,7.8); type2i4(&x,8) -> 3
223{
224 r_4* x4; r_8* x8; int_4 lw=0;
225 if(nbytes==4) {x4 = (r_4*)val; lw = (int_4)(*x4);}
226 else {x8 = (r_8*)val; lw = (int_4)(*x8);}
227 return lw;
228}
229
230////////////////////////////////////////////////////////////////////////////////////
231//! Interface to Lapack linear system solver driver s/d/c/zgesv().
232/*! Solve the linear system a * x = b using LU factorization.
233 Input arrays should have FortranMemory mapping (column packed).
234 \param a : input matrix, overwritten on output
235 \param b : input-output, input vector b, contains x on exit
236 \return : return code from lapack driver _gesv()
237 */
238template <class T>
239int LapackServer<T>::LinSolve(TArray<T>& a, TArray<T> & b)
240{
241 if ( ( a.NbDimensions() != 2 ) || ( b.NbDimensions() != 2 ) )
242 throw(SzMismatchError("LapackServer::LinSolve(a,b) a Or b NbDimensions() != 2"));
243
244 int_4 rowa = a.RowsKA();
245 int_4 cola = a.ColsKA();
246 int_4 rowb = b.RowsKA();
247 int_4 colb = b.ColsKA();
248 if ( a.Size(rowa) != a.Size(cola))
249 throw(SzMismatchError("LapackServer::LinSolve(a,b) a Not a square Array"));
250 if ( a.Size(rowa) != b.Size(rowb))
251 throw(SzMismatchError("LapackServer::LinSolve(a,b) RowSize(a <> b) "));
252
253 if (!a.IsPacked(rowa) || !b.IsPacked(rowb))
254 throw(SzMismatchError("LapackServer::LinSolve(a,b) a Or b Not Column Packed"));
255
256 int_4 n = a.Size(rowa);
257 int_4 nrhs = b.Size(colb);
258 int_4 lda = a.Step(cola);
259 int_4 ldb = b.Step(colb);
260 int_4 info;
261 int_4* ipiv = new int_4[n];
262
263 if (typeid(T) == typeid(r_4) )
264 sgesv_(&n, &nrhs, (r_4 *)a.Data(), &lda, ipiv, (r_4 *)b.Data(), &ldb, &info);
265 else if (typeid(T) == typeid(r_8) )
266 dgesv_(&n, &nrhs, (r_8 *)a.Data(), &lda, ipiv, (r_8 *)b.Data(), &ldb, &info);
267 else if (typeid(T) == typeid(complex<r_4>) )
268 cgesv_(&n, &nrhs, (complex<r_4> *)a.Data(), &lda, ipiv,
269 (complex<r_4> *)b.Data(), &ldb, &info);
270 else if (typeid(T) == typeid(complex<r_8>) )
271 zgesv_(&n, &nrhs, (complex<r_8> *)a.Data(), &lda, ipiv,
272 (complex<r_8> *)b.Data(), &ldb, &info);
273 else {
274 delete[] ipiv;
275 string tn = typeid(T).name();
276 cerr << " LapackServer::LinSolve(a,b) - Unsupported DataType T = " << tn << endl;
277 throw TypeMismatchExc("LapackServer::LinSolve(a,b) - Unsupported DataType (T)");
278 }
279 delete[] ipiv;
280 return(info);
281}
282
283////////////////////////////////////////////////////////////////////////////////////
284//! Interface to Lapack linear system solver driver s/d/c/zsysv().
285/*! Solve the linear system a * x = b with a symetric matrix using LU factorization.
286 Input arrays should have FortranMemory mapping (column packed).
287 \param a : input matrix symetric , overwritten on output
288 \param b : input-output, input vector b, contains x on exit
289 \return : return code from lapack driver _sysv()
290 */
291template <class T>
292int LapackServer<T>::LinSolveSym(TArray<T>& a, TArray<T> & b)
293// --- REMARQUES DE CMV ---
294// 1./ contrairement a ce qui est dit dans la doc, il s'agit
295// de matrices SYMETRIQUES complexes et non HERMITIENNES !!!
296// 2./ pourquoi les routines de LinSolve pour des matrices symetriques
297// sont plus de deux fois plus lentes que les LinSolve generales sur OSF
298// et sensiblement plus lentes sous Linux ???
299{
300 if ( ( a.NbDimensions() != 2 ) || ( b.NbDimensions() != 2 ) )
301 throw(SzMismatchError("LapackServer::LinSolveSym(a,b) a Or b NbDimensions() != 2"));
302 int_4 rowa = a.RowsKA();
303 int_4 cola = a.ColsKA();
304 int_4 rowb = b.RowsKA();
305 int_4 colb = b.ColsKA();
306 if ( a.Size(rowa) != a.Size(cola))
307 throw(SzMismatchError("LapackServer::LinSolveSym(a,b) a Not a square Array"));
308 if ( a.Size(rowa) != b.Size(rowb))
309 throw(SzMismatchError("LapackServer::LinSolveSym(a,b) RowSize(a <> b) "));
310
311 if (!a.IsPacked(rowa) || !b.IsPacked(rowb))
312 throw(SzMismatchError("LapackServer::LinSolveSym(a,b) a Or b Not Column Packed"));
313
314 int_4 n = a.Size(rowa);
315 int_4 nrhs = b.Size(colb);
316 int_4 lda = a.Step(cola);
317 int_4 ldb = b.Step(colb);
318 int_4 info = 0;
319 int_4* ipiv = new int_4[n];
320 int_4 lwork = -1;
321 T * work = NULL;
322 T wkget[2];
323
324 char uplo = 'U'; // char uplo = 'L';
325 char struplo[5]; struplo[0] = uplo; struplo[1] = '\0';
326
327 if (typeid(T) == typeid(r_4) ) {
328 ssysv_(&uplo, &n, &nrhs, (r_4 *)a.Data(), &lda, ipiv, (r_4 *)b.Data(), &ldb,
329 (r_4 *)wkget, &lwork, &info);
330 lwork = type2i4(&wkget[0],4); work = new T[lwork +GARDMEM];
331 ssysv_(&uplo, &n, &nrhs, (r_4 *)a.Data(), &lda, ipiv, (r_4 *)b.Data(), &ldb,
332 (r_4 *)work, &lwork, &info);
333 } else if (typeid(T) == typeid(r_8) ) {
334 dsysv_(&uplo, &n, &nrhs, (r_8 *)a.Data(), &lda, ipiv, (r_8 *)b.Data(), &ldb,
335 (r_8 *)wkget, &lwork, &info);
336 lwork = type2i4(&wkget[0],8); work = new T[lwork +GARDMEM];
337 dsysv_(&uplo, &n, &nrhs, (r_8 *)a.Data(), &lda, ipiv, (r_8 *)b.Data(), &ldb,
338 (r_8 *)work, &lwork, &info);
339 } else if (typeid(T) == typeid(complex<r_4>) ) {
340 csysv_(&uplo, &n, &nrhs, (complex<r_4> *)a.Data(), &lda, ipiv,
341 (complex<r_4> *)b.Data(), &ldb,
342 (complex<r_4> *)wkget, &lwork, &info);
343 lwork = type2i4(&wkget[0],4); work = new T[lwork +GARDMEM];
344 csysv_(&uplo, &n, &nrhs, (complex<r_4> *)a.Data(), &lda, ipiv,
345 (complex<r_4> *)b.Data(), &ldb,
346 (complex<r_4> *)work, &lwork, &info);
347 } else if (typeid(T) == typeid(complex<r_8>) ) {
348 zsysv_(&uplo, &n, &nrhs, (complex<r_8> *)a.Data(), &lda, ipiv,
349 (complex<r_8> *)b.Data(), &ldb,
350 (complex<r_8> *)wkget, &lwork, &info);
351 lwork = type2i4(&wkget[0],8); work = new T[lwork +GARDMEM];
352 zsysv_(&uplo, &n, &nrhs, (complex<r_8> *)a.Data(), &lda, ipiv,
353 (complex<r_8> *)b.Data(), &ldb,
354 (complex<r_8> *)work, &lwork, &info);
355 } else {
356 if(work) delete[] work;
357 delete[] ipiv;
358 string tn = typeid(T).name();
359 cerr << " LapackServer::LinSolveSym(a,b) - Unsupported DataType T = " << tn << endl;
360 throw TypeMismatchExc("LapackServer::LinSolveSym(a,b) - Unsupported DataType (T)");
361 }
362 if(work) delete[] work;
363 delete[] ipiv;
364 return(info);
365}
366
367////////////////////////////////////////////////////////////////////////////////////
368//! Interface to Lapack least squares solver driver s/d/c/zgels().
369/*! Solves the linear least squares problem defined by an m-by-n matrix
370 \b a and an m element vector \b b , using QR or LQ factorization .
371 A solution \b x to the overdetermined system of linear equations
372 b = a * x is computed, minimizing the norm of b-a*x.
373 Underdetermined systems (m<n) are not yet handled.
374 Inout arrays should have FortranMemory mapping (column packed).
375 \param a : input matrix, overwritten on output
376 \param b : input-output, input vector b, contains x on exit.
377 \return : return code from lapack driver _gels()
378 \warning : b is not resized.
379 */
380/*
381 $CHECK$ - A faire - cas m<n
382 If the linear system is underdetermined, the minimum norm
383 solution is computed.
384*/
385
386template <class T>
387int LapackServer<T>::LeastSquareSolve(TArray<T>& a, TArray<T> & b)
388{
389 if ( ( a.NbDimensions() != 2 ) || ( b.NbDimensions() != 2 ) )
390 throw(SzMismatchError("LapackServer::LeastSquareSolve(a,b) a Or b NbDimensions() != 2"));
391
392 int_4 rowa = a.RowsKA();
393 int_4 cola = a.ColsKA();
394 int_4 rowb = b.RowsKA();
395 int_4 colb = b.ColsKA();
396
397
398 if ( a.Size(rowa) != b.Size(rowb))
399 throw(SzMismatchError("LapackServer::LeastSquareSolve(a,b) RowSize(a <> b) "));
400
401 if (!a.IsPacked(rowa) || !b.IsPacked(rowb))
402 throw(SzMismatchError("LapackServer::LeastSquareSolve(a,b) a Or b Not Column Packed"));
403
404 if ( a.Size(rowa) < a.Size(cola)) { // $CHECK$ - m<n a changer
405 cout << " LapackServer<T>::LeastSquareSolve() - m<n - Not yet implemented for "
406 << " underdetermined systems ! " << endl;
407 throw(SzMismatchError("LapackServer::LeastSquareSolve(a,b) NRows<NCols - "));
408 }
409 int_4 m = a.Size(rowa);
410 int_4 n = a.Size(cola);
411 int_4 nrhs = b.Size(colb);
412
413 int_4 lda = a.Step(cola);
414 int_4 ldb = b.Step(colb);
415 int_4 info;
416
417 int_4 minmn = (m < n) ? m : n;
418 int_4 maxmn = (m > n) ? m : n;
419 int_4 maxmnrhs = (nrhs > maxmn) ? nrhs : maxmn;
420 if (maxmnrhs < 1) maxmnrhs = 1;
421
422 int_4 lwork = -1; //minmn+maxmnrhs*5;
423 T * work = NULL;
424 T wkget[2];
425
426 char trans = 'N';
427
428 if (typeid(T) == typeid(r_4) ) {
429 sgels_(&trans, &m, &n, &nrhs, (r_4 *)a.Data(), &lda,
430 (r_4 *)b.Data(), &ldb, (r_4 *)wkget, &lwork, &info);
431 lwork = type2i4(&wkget[0],4); work = new T[lwork +GARDMEM];
432 sgels_(&trans, &m, &n, &nrhs, (r_4 *)a.Data(), &lda,
433 (r_4 *)b.Data(), &ldb, (r_4 *)work, &lwork, &info);
434 } else if (typeid(T) == typeid(r_8) ) {
435 dgels_(&trans, &m, &n, &nrhs, (r_8 *)a.Data(), &lda,
436 (r_8 *)b.Data(), &ldb, (r_8 *)wkget, &lwork, &info);
437 lwork = type2i4(&wkget[0],8); work = new T[lwork +GARDMEM];
438 dgels_(&trans, &m, &n, &nrhs, (r_8 *)a.Data(), &lda,
439 (r_8 *)b.Data(), &ldb, (r_8 *)work, &lwork, &info);
440 } else if (typeid(T) == typeid(complex<r_4>) ) {
441 cgels_(&trans, &m, &n, &nrhs, (complex<r_4> *)a.Data(), &lda,
442 (complex<r_4> *)b.Data(), &ldb, (complex<r_4> *)wkget, &lwork, &info);
443 lwork = type2i4(&wkget[0],4); work = new T[lwork +GARDMEM];
444 cgels_(&trans, &m, &n, &nrhs, (complex<r_4> *)a.Data(), &lda,
445 (complex<r_4> *)b.Data(), &ldb, (complex<r_4> *)work, &lwork, &info);
446 } else if (typeid(T) == typeid(complex<r_8>) ) {
447 zgels_(&trans, &m, &n, &nrhs, (complex<r_8> *)a.Data(), &lda,
448 (complex<r_8> *)b.Data(), &ldb, (complex<r_8> *)wkget, &lwork, &info);
449 lwork = type2i4(&wkget[0],8); work = new T[lwork +GARDMEM];
450 zgels_(&trans, &m, &n, &nrhs, (complex<r_8> *)a.Data(), &lda,
451 (complex<r_8> *)b.Data(), &ldb, (complex<r_8> *)work, &lwork, &info);
452 } else {
453 if(work) delete [] work; work=NULL;
454 string tn = typeid(T).name();
455 cerr << " LapackServer::LeastSquareSolve(a,b) - Unsupported DataType T = " << tn << endl;
456 throw TypeMismatchExc("LapackServer::LeastSquareSolve(a,b) - Unsupported DataType (T)");
457 }
458 if(work) delete [] work;
459 return(info);
460}
461
462////////////////////////////////////////////////////////////////////////////////////
463//! Square matrix inversion using Lapack linear system solver
464/*! Compute the inverse of a square matrix using linear system solver routine
465 Input arrays should have FortranMemory mapping (column packed).
466 \param a : input matrix, overwritten on output
467 \param ainv : output matrix, contains inverse(a) on exit.
468 ainv is allocated if it has size 0
469 If not allocated, ainv is automatically
470 \return : return code from LapackServer::LinSolve()
471 \sa LapackServer::LinSolve()
472 */
473template <class T>
474int LapackServer<T>::ComputeInverse(TMatrix<T>& a, TMatrix<T> & ainv)
475{
476 if ( a.NbDimensions() != 2 )
477 throw(SzMismatchError("LapackServer::Inverse() NDim(a) != 2"));
478 if ( a.GetMemoryMapping() != BaseArray::FortranMemoryMapping )
479 throw(SzMismatchError("LapackServer::Inverse() a NOT in FortranMemoryMapping"));
480 if ( a.NRows() != a.NCols() )
481 throw(SzMismatchError("LapackServer::Inverse() a NOT square matrix (a.NRows!=a.NCols)"));
482 if (ainv.IsAllocated()) {
483 bool smo, ssz;
484 ssz = a.CompareSizes(ainv, smo);
485 if ( (ssz == false) || (smo == false) )
486 throw(SzMismatchError("LapackServer::Inverse() ainv<>a Size/MemOrg mismatch "));
487 }
488 else ainv.SetSize(a.NRows(), a.NCols(), BaseArray::FortranMemoryMapping, false);
489 ainv = IdentityMatrix();
490 return LinSolve(a, ainv);
491}
492
493////////////////////////////////////////////////////////////////////////////////////
494//! Interface to Lapack least squares solver driver s/d/c/zgelsd().
495/*! Solves the linear least squares problem defined by an m-by-n matrix
496 \b a and an m element vector \b b , using SVD factorization Divide and Conquer.
497 Inout arrays should have FortranMemory mapping (column packed).
498 \param rcond : definition of zero value (S(i) <= RCOND*S(0) are treated as zero).
499 If RCOND < 0, machine precision is used instead.
500 \param a : input matrix, overwritten on output
501 \param b : input vector b overwritten by solution on output (beware of size changing)
502 \param x : output matrix of solutions.
503 \param rank : output the rank of the matrix.
504 \return : return code from lapack driver _gelsd()
505 \warning : b is not resized.
506 */
507template <class T>
508int LapackServer<T>::LeastSquareSolveSVD_DC(TMatrix<T>& a,TMatrix<T>& b,TVector<r_8>& s,int_4& rank,r_8 rcond)
509{
510 if ( ( a.NbDimensions() != 2 ) )
511 throw(SzMismatchError("LapackServer::LeastSquareSolveSVD_DC(a,b) a != 2"));
512
513 if (!a.IsPacked() || !b.IsPacked())
514 throw(SzMismatchError("LapackServer::LeastSquareSolveSVD_DC(a,b) a Or b Not Packed"));
515
516 int_4 m = a.NRows();
517 int_4 n = a.NCols();
518
519 if(b.NRows() != m)
520 throw(SzMismatchError("LapackServer::LeastSquareSolveSVD_DC(a,b) bad matching dim between a and b"));
521
522 int_4 nrhs = b.NCols();
523 int_4 minmn = (m < n) ? m : n;
524 int_4 maxmn = (m > n) ? m : n;
525
526 int_4 lda = m;
527 int_4 ldb = maxmn;
528 int_4 info;
529
530 { // Use {} for automatic des-allocation of "bsave"
531 TMatrix<T> bsave(m,nrhs); bsave.SetMemoryMapping(BaseArray::FortranMemoryMapping);
532 bsave = b;
533 b.ReSize(maxmn,nrhs); b = (T) 0;
534 for(int i=0;i<m;i++) for(int j=0;j<nrhs;j++) b(i,j) = bsave(i,j);
535 } // Use {} for automatic des-allocation of "bsave"
536 s.ReSize(minmn);
537
538 int_4 smlsiz = 25; // Normallement ilaenv_en_C(9,...) renvoie toujours 25
539 if(typeid(T) == typeid(r_4) ) smlsiz = ilaenv_en_C(9,"SGELSD"," ",0,0,0,0);
540 else if(typeid(T) == typeid(r_8) ) smlsiz = ilaenv_en_C(9,"DGELSD"," ",0,0,0,0);
541 else if(typeid(T) == typeid(complex<r_4>) ) smlsiz = ilaenv_en_C(9,"CGELSD"," ",0,0,0,0);
542 else if(typeid(T) == typeid(complex<r_8>) ) smlsiz = ilaenv_en_C(9,"ZGELSD"," ",0,0,0,0);
543 if(smlsiz<0) smlsiz = 0;
544 r_8 dum = log((r_8)minmn/(r_8)(smlsiz+1.)) / log(2.);
545 int_4 nlvl = int_4(dum) + 1; if(nlvl<0) nlvl = 0;
546
547 T * work = NULL;
548 int_4 * iwork = NULL;
549 int_4 lwork=-1, lrwork;
550 T wkget[2];
551
552 if(typeid(T) == typeid(r_4) ) {
553 r_4* sloc = new r_4[minmn];
554 r_4 srcond = rcond;
555 iwork = new int_4[3*minmn*nlvl+11*minmn +GARDMEM];
556 sgelsd_(&m,&n,&nrhs,(r_4*)a.Data(),&lda,
557 (r_4*)b.Data(),&ldb,(r_4*)sloc,&srcond,&rank,
558 (r_4*)wkget,&lwork,(int_4*)iwork,&info);
559 lwork = type2i4(&wkget[0],4); work = new T[lwork +GARDMEM];
560 sgelsd_(&m,&n,&nrhs,(r_4*)a.Data(),&lda,
561 (r_4*)b.Data(),&ldb,(r_4*)sloc,&srcond,&rank,
562 (r_4*)work,&lwork,(int_4*)iwork,&info);
563 for(int_4 i=0;i<minmn;i++) s(i) = sloc[i];
564 delete [] sloc;
565 } else if(typeid(T) == typeid(r_8) ) {
566 iwork = new int_4[3*minmn*nlvl+11*minmn +GARDMEM];
567 dgelsd_(&m,&n,&nrhs,(r_8*)a.Data(),&lda,
568 (r_8*)b.Data(),&ldb,(r_8*)s.Data(),&rcond,&rank,
569 (r_8*)wkget,&lwork,(int_4*)iwork,&info);
570 lwork = type2i4(&wkget[0],8); work = new T[lwork +GARDMEM];
571 dgelsd_(&m,&n,&nrhs,(r_8*)a.Data(),&lda,
572 (r_8*)b.Data(),&ldb,(r_8*)s.Data(),&rcond,&rank,
573 (r_8*)work,&lwork,(int_4*)iwork,&info);
574 } else if(typeid(T) == typeid(complex<r_4>) ) {
575 // Cf meme remarque que ci-dessous (complex<r_8)
576 lrwork = 10*minmn + 2*minmn*smlsiz + 8*minmn*nlvl + 3*smlsiz*nrhs + (smlsiz+1)*(smlsiz+1);
577 int_4 lrwork_d = 12*minmn + 2*minmn*smlsiz + 8*minmn*nlvl + minmn*nrhs + (smlsiz+1)*(smlsiz+1);
578 if(lrwork_d > lrwork) lrwork = lrwork_d;
579 r_4* rwork = new r_4[lrwork +GARDMEM];
580 iwork = new int_4[3*minmn*nlvl+11*minmn +GARDMEM];
581 r_4* sloc = new r_4[minmn];
582 r_4 srcond = rcond;
583 cgelsd_(&m,&n,&nrhs,(complex<r_4>*)a.Data(),&lda,
584 (complex<r_4>*)b.Data(),&ldb,(r_4*)sloc,&srcond,&rank,
585 (complex<r_4>*)wkget,&lwork,(r_4*)rwork,(int_4*)iwork,&info);
586 lwork = type2i4(&wkget[0],4); work = new T[lwork +GARDMEM];
587 cgelsd_(&m,&n,&nrhs,(complex<r_4>*)a.Data(),&lda,
588 (complex<r_4>*)b.Data(),&ldb,(r_4*)sloc,&srcond,&rank,
589 (complex<r_4>*)work,&lwork,(r_4*)rwork,(int_4*)iwork,&info);
590 for(int_4 i=0;i<minmn;i++) s(i) = sloc[i];
591 delete [] sloc; delete [] rwork;
592 } else if(typeid(T) == typeid(complex<r_8>) ) {
593 // CMV: Bizarrement, la formule donnee dans zgelsd() plante pour des N grands (500)
594 // On prend (par analogie) la formule pour "lwork" de dgelsd()
595 lrwork = 10*minmn + 2*minmn*smlsiz + 8*minmn*nlvl + 3*smlsiz*nrhs + (smlsiz+1)*(smlsiz+1);
596 int_4 lrwork_d = 12*minmn + 2*minmn*smlsiz + 8*minmn*nlvl + minmn*nrhs + (smlsiz+1)*(smlsiz+1);
597 if(lrwork_d > lrwork) lrwork = lrwork_d;
598 r_8* rwork = new r_8[lrwork +GARDMEM];
599 iwork = new int_4[3*minmn*nlvl+11*minmn +GARDMEM];
600 zgelsd_(&m,&n,&nrhs,(complex<r_8>*)a.Data(),&lda,
601 (complex<r_8>*)b.Data(),&ldb,(r_8*)s.Data(),&rcond,&rank,
602 (complex<r_8>*)wkget,&lwork,(r_8*)rwork,(int_4*)iwork,&info);
603 lwork = type2i4(&wkget[0],8); work = new T[lwork +GARDMEM];
604 zgelsd_(&m,&n,&nrhs,(complex<r_8>*)a.Data(),&lda,
605 (complex<r_8>*)b.Data(),&ldb,(r_8*)s.Data(),&rcond,&rank,
606 (complex<r_8>*)work,&lwork,(r_8*)rwork,(int_4*)iwork,&info);
607 delete [] rwork;
608 } else {
609 if(work) delete [] work; work=NULL;
610 if(iwork) delete [] iwork; iwork=NULL;
611 string tn = typeid(T).name();
612 cerr << " LapackServer::LeastSquareSolveSVD_DC(a,b) - Unsupported DataType T = " << tn << endl;
613 throw TypeMismatchExc("LapackServer::LeastSquareSolveSVD_DC(a,b) - Unsupported DataType (T)");
614 }
615
616 if(work) delete [] work; if(iwork) delete [] iwork;
617 return(info);
618}
619
620
621////////////////////////////////////////////////////////////////////////////////////
622//! Interface to Lapack SVD driver s/d/c/zgesv().
623/*! Computes the vector of singular values of \b a. Input arrays
624 should have FortranMemoryMapping (column packed).
625 \param a : input m-by-n matrix
626 \param s : Vector of min(m,n) singular values (descending order)
627 \return : return code from lapack driver _gesvd()
628 */
629
630template <class T>
631int LapackServer<T>::SVD(TArray<T>& a, TArray<T> & s)
632{
633 return (SVDDriver(a, s, NULL, NULL) );
634}
635
636//! Interface to Lapack SVD driver s/d/c/zgesv().
637/*! Computes the vector of singular values of \b a, as well as
638 right and left singular vectors of \b a.
639 \f[
640 A = U \Sigma V^T , ( A = U \Sigma V^H \ complex)
641 \f]
642 \f[
643 A v_i = \sigma_i u_i \ and A^T u_i = \sigma_i v_i \ (A^H \ complex)
644 \f]
645 U and V are orthogonal (unitary) matrices.
646 \param a : input m-by-n matrix (in FortranMemoryMapping)
647 \param s : Vector of min(m,n) singular values (descending order)
648 \param u : m-by-m Matrix of left singular vectors
649 \param vt : Transpose of right singular vectors (n-by-n matrix).
650 \return : return code from lapack driver _gesvd()
651 */
652template <class T>
653int LapackServer<T>::SVD(TArray<T>& a, TArray<T> & s, TArray<T> & u, TArray<T> & vt)
654{
655 return (SVDDriver(a, s, &u, &vt) );
656}
657
658
659//! Interface to Lapack SVD driver s/d/c/zgesv().
660template <class T>
661int LapackServer<T>::SVDDriver(TArray<T>& a, TArray<T> & s, TArray<T>* up, TArray<T>* vtp)
662{
663 if ( ( a.NbDimensions() != 2 ) )
664 throw(SzMismatchError("LapackServer::SVDDriver(a, ...) a.NbDimensions() != 2"));
665
666 int_4 rowa = a.RowsKA();
667 int_4 cola = a.ColsKA();
668
669 if ( !a.IsPacked(rowa) )
670 throw(SzMismatchError("LapackServer::SVDDriver(a, ...) a Not Column Packed "));
671
672 int_4 m = a.Size(rowa);
673 int_4 n = a.Size(cola);
674 int_4 maxmn = (m > n) ? m : n;
675 int_4 minmn = (m < n) ? m : n;
676
677 char jobu, jobvt;
678 jobu = 'N';
679 jobvt = 'N';
680
681 sa_size_t sz[2];
682 if ( up != NULL) {
683 if ( dynamic_cast< TVector<T> * > (vtp) )
684 throw( TypeMismatchExc("LapackServer::SVDDriver() Wrong type (=TVector<T>) for u !") );
685 up->SetMemoryMapping(BaseArray::FortranMemoryMapping);
686 sz[0] = sz[1] = m;
687 up->ReSize(2, sz );
688 jobu = 'A';
689 }
690 else {
691 up = new TMatrix<T>(1,1);
692 jobu = 'N';
693 }
694 if ( vtp != NULL) {
695 if ( dynamic_cast< TVector<T> * > (vtp) )
696 throw( TypeMismatchExc("LapackServer::SVDDriver() Wrong type (=TVector<T>) for vt !") );
697 vtp->SetMemoryMapping(BaseArray::FortranMemoryMapping);
698 sz[0] = sz[1] = n;
699 vtp->ReSize(2, sz );
700 jobvt = 'A';
701 }
702 else {
703 vtp = new TMatrix<T>(1,1);
704 jobvt = 'N';
705 }
706
707 TVector<T> *vs = dynamic_cast< TVector<T> * > (&s);
708 if (vs) vs->ReSize(minmn);
709 else {
710 TMatrix<T> *ms = dynamic_cast< TMatrix<T> * > (&s);
711 if (ms) ms->ReSize(minmn,1);
712 else {
713 sz[0] = minmn; sz[1] = 1;
714 s.ReSize(1, sz);
715 }
716 }
717
718 int_4 lda = a.Step(a.ColsKA());
719 int_4 ldu = up->Step(up->ColsKA());
720 int_4 ldvt = vtp->Step(vtp->ColsKA());
721 int_4 info;
722
723 int_4 lwork = -1; // maxmn*5 *wspace_size_factor;
724 T * work = NULL; // = new T[lwork];
725 T wkget[2];
726
727 if (typeid(T) == typeid(r_4) ) {
728 sgesvd_(&jobu, &jobvt, &m, &n, (r_4 *)a.Data(), &lda,
729 (r_4 *)s.Data(), (r_4 *) up->Data(), &ldu, (r_4 *)vtp->Data(), &ldvt,
730 (r_4 *)wkget, &lwork, &info);
731 lwork = type2i4(&wkget[0],4); work = new T[lwork +GARDMEM];
732 sgesvd_(&jobu, &jobvt, &m, &n, (r_4 *)a.Data(), &lda,
733 (r_4 *)s.Data(), (r_4 *) up->Data(), &ldu, (r_4 *)vtp->Data(), &ldvt,
734 (r_4 *)work, &lwork, &info);
735 } else if (typeid(T) == typeid(r_8) ) {
736 dgesvd_(&jobu, &jobvt, &m, &n, (r_8 *)a.Data(), &lda,
737 (r_8 *)s.Data(), (r_8 *) up->Data(), &ldu, (r_8 *)vtp->Data(), &ldvt,
738 (r_8 *)wkget, &lwork, &info);
739 lwork = type2i4(&wkget[0],8); work = new T[lwork +GARDMEM];
740 dgesvd_(&jobu, &jobvt, &m, &n, (r_8 *)a.Data(), &lda,
741 (r_8 *)s.Data(), (r_8 *) up->Data(), &ldu, (r_8 *)vtp->Data(), &ldvt,
742 (r_8 *)work, &lwork, &info);
743 } else if (typeid(T) == typeid(complex<r_4>) ) {
744 r_4 * rwork = new r_4[5*minmn +GARDMEM];
745 r_4 * sloc = new r_4[minmn];
746 cgesvd_(&jobu, &jobvt, &m, &n, (complex<r_4> *)a.Data(), &lda,
747 (r_4 *)sloc, (complex<r_4> *) up->Data(), &ldu,
748 (complex<r_4> *)vtp->Data(), &ldvt,
749 (complex<r_4> *)wkget, &lwork, (r_4 *)rwork, &info);
750 lwork = type2i4(&wkget[0],4); work = new T[lwork +GARDMEM];
751 cgesvd_(&jobu, &jobvt, &m, &n, (complex<r_4> *)a.Data(), &lda,
752 (r_4 *)sloc, (complex<r_4> *) up->Data(), &ldu,
753 (complex<r_4> *)vtp->Data(), &ldvt,
754 (complex<r_4> *)work, &lwork, (r_4 *)rwork, &info);
755 for(int_4 i=0;i<minmn;i++) s[i] = sloc[i];
756 delete [] rwork; delete [] sloc;
757 } else if (typeid(T) == typeid(complex<r_8>) ) {
758 r_8 * rwork = new r_8[5*minmn +GARDMEM];
759 r_8 * sloc = new r_8[minmn];
760 zgesvd_(&jobu, &jobvt, &m, &n, (complex<r_8> *)a.Data(), &lda,
761 (r_8 *)sloc, (complex<r_8> *) up->Data(), &ldu,
762 (complex<r_8> *)vtp->Data(), &ldvt,
763 (complex<r_8> *)wkget, &lwork, (r_8 *)rwork, &info);
764 lwork = type2i4(&wkget[0],8); work = new T[lwork +GARDMEM];
765 zgesvd_(&jobu, &jobvt, &m, &n, (complex<r_8> *)a.Data(), &lda,
766 (r_8 *)sloc, (complex<r_8> *) up->Data(), &ldu,
767 (complex<r_8> *)vtp->Data(), &ldvt,
768 (complex<r_8> *)work, &lwork, (r_8 *)rwork, &info);
769 for(int_4 i=0;i<minmn;i++) s[i] = sloc[i];
770 delete [] rwork; delete [] sloc;
771 } else {
772 if(work) delete [] work; work=NULL;
773 if (jobu == 'N') delete up;
774 if (jobvt == 'N') delete vtp;
775 string tn = typeid(T).name();
776 cerr << " LapackServer::SVDDriver(...) - Unsupported DataType T = " << tn << endl;
777 throw TypeMismatchExc("LapackServer::SVDDriver(a,b) - Unsupported DataType (T)");
778 }
779
780 if(work) delete [] work;
781 if (jobu == 'N') delete up;
782 if (jobvt == 'N') delete vtp;
783 return(info);
784}
785
786
787//! Interface to Lapack SVD driver s/d/c/zgesdd().
788/*! Same as SVD but with Divide and Conquer method */
789template <class T>
790int LapackServer<T>::SVD_DC(TMatrix<T>& a, TVector<r_8>& s, TMatrix<T>& u, TMatrix<T>& vt)
791{
792
793 if ( !a.IsPacked() )
794 throw(SzMismatchError("LapackServer::SVD_DC(a, ...) a Not Packed "));
795
796 int_4 m = a.NRows();
797 int_4 n = a.NCols();
798 int_4 maxmn = (m > n) ? m : n;
799 int_4 minmn = (m < n) ? m : n;
800 int_4 supermax = 4*minmn*minmn+4*minmn; if(maxmn>supermax) supermax=maxmn;
801
802 char jobz = 'A';
803
804 s.ReSize(minmn);
805 u.ReSize(m,m);
806 vt.ReSize(n,n);
807
808 int_4 lda = m;
809 int_4 ldu = m;
810 int_4 ldvt = n;
811 int_4 info;
812 int_4 lwork=-1;
813 T * work = NULL;
814 int_4 * iwork = NULL;
815 T wkget[2];
816
817 if(typeid(T) == typeid(r_4) ) {
818 r_4* sloc = new r_4[minmn];
819 iwork = new int_4[8*minmn +GARDMEM];
820 sgesdd_(&jobz,&m,&n,(r_4*)a.Data(),&lda,
821 (r_4*)sloc,(r_4*)u.Data(),&ldu,(r_4*)vt.Data(),&ldvt,
822 (r_4*)wkget,&lwork,(int_4*)iwork,&info);
823 lwork = type2i4(&wkget[0],4); work = new T[lwork +GARDMEM];
824 sgesdd_(&jobz,&m,&n,(r_4*)a.Data(),&lda,
825 (r_4*)sloc,(r_4*)u.Data(),&ldu,(r_4*)vt.Data(),&ldvt,
826 (r_4*)work,&lwork,(int_4*)iwork,&info);
827 for(int_4 i=0;i<minmn;i++) s(i) = (r_8) sloc[i];
828 delete [] sloc;
829 } else if(typeid(T) == typeid(r_8) ) {
830 iwork = new int_4[8*minmn +GARDMEM];
831 dgesdd_(&jobz,&m,&n,(r_8*)a.Data(),&lda,
832 (r_8*)s.Data(),(r_8*)u.Data(),&ldu,(r_8*)vt.Data(),&ldvt,
833 (r_8*)wkget,&lwork,(int_4*)iwork,&info);
834 lwork = type2i4(&wkget[0],8); work = new T[lwork +GARDMEM];
835 dgesdd_(&jobz,&m,&n,(r_8*)a.Data(),&lda,
836 (r_8*)s.Data(),(r_8*)u.Data(),&ldu,(r_8*)vt.Data(),&ldvt,
837 (r_8*)work,&lwork,(int_4*)iwork,&info);
838 } else if(typeid(T) == typeid(complex<r_4>) ) {
839 r_4* sloc = new r_4[minmn];
840 r_4* rwork = new r_4[5*minmn*minmn+5*minmn +GARDMEM];
841 iwork = new int_4[8*minmn +GARDMEM];
842 cgesdd_(&jobz,&m,&n,(complex<r_4>*)a.Data(),&lda,
843 (r_4*)sloc,(complex<r_4>*)u.Data(),&ldu,(complex<r_4>*)vt.Data(),&ldvt,
844 (complex<r_4>*)wkget,&lwork,(r_4*)rwork,(int_4*)iwork,&info);
845 lwork = type2i4(&wkget[0],4); work = new T[lwork +GARDMEM];
846 cgesdd_(&jobz,&m,&n,(complex<r_4>*)a.Data(),&lda,
847 (r_4*)sloc,(complex<r_4>*)u.Data(),&ldu,(complex<r_4>*)vt.Data(),&ldvt,
848 (complex<r_4>*)work,&lwork,(r_4*)rwork,(int_4*)iwork,&info);
849 for(int_4 i=0;i<minmn;i++) s(i) = (r_8) sloc[i];
850 delete [] sloc; delete [] rwork;
851 } else if(typeid(T) == typeid(complex<r_8>) ) {
852 r_8* rwork = new r_8[5*minmn*minmn+5*minmn +GARDMEM];
853 iwork = new int_4[8*minmn +GARDMEM];
854 zgesdd_(&jobz,&m,&n,(complex<r_8>*)a.Data(),&lda,
855 (r_8*)s.Data(),(complex<r_8>*)u.Data(),&ldu,(complex<r_8>*)vt.Data(),&ldvt,
856 (complex<r_8>*)wkget,&lwork,(r_8*)rwork,(int_4*)iwork,&info);
857 lwork = type2i4(&wkget[0],8); work = new T[lwork +GARDMEM];
858 zgesdd_(&jobz,&m,&n,(complex<r_8>*)a.Data(),&lda,
859 (r_8*)s.Data(),(complex<r_8>*)u.Data(),&ldu,(complex<r_8>*)vt.Data(),&ldvt,
860 (complex<r_8>*)work,&lwork,(r_8*)rwork,(int_4*)iwork,&info);
861 delete [] rwork;
862 } else {
863 if(work) delete [] work; work=NULL;
864 if(iwork) delete [] iwork; iwork=NULL;
865 string tn = typeid(T).name();
866 cerr << " LapackServer::SVD_DC(...) - Unsupported DataType T = " << tn << endl;
867 throw TypeMismatchExc("LapackServer::SVD_DC - Unsupported DataType (T)");
868 }
869
870 if(work) delete [] work; if(iwork) delete [] iwork;
871 return(info);
872}
873
874
875////////////////////////////////////////////////////////////////////////////////////
876/*! Computes the eigen values and eigen vectors of a symetric (or hermitian) matrix \b a.
877 Input arrays should have FortranMemoryMapping (column packed).
878 \param a : input symetric (or hermitian) n-by-n matrix
879 \param b : Vector of eigenvalues (descending order)
880 \param eigenvector : if true compute eigenvectors, if not only eigenvalues
881 \param a : on return array of eigenvectors (same order than eval, one vector = one column)
882 \return : return code from lapack driver
883 */
884
885template <class T>
886int LapackServer<T>::LapackEigenSym(TArray<T>& a, TVector<r_8>& b, bool eigenvector)
887{
888 if ( a.NbDimensions() != 2 )
889 throw(SzMismatchError("LapackServer::LapackEigenSym(a,b) a NbDimensions() != 2"));
890 int_4 rowa = a.RowsKA();
891 int_4 cola = a.ColsKA();
892 if ( a.Size(rowa) != a.Size(cola))
893 throw(SzMismatchError("LapackServer::LapackEigenSym(a,b) a Not a square Array"));
894 if (!a.IsPacked(rowa))
895 throw(SzMismatchError("LapackServer::LapackEigenSym(a,b) a Not Column Packed"));
896
897 char uplo='U';
898 char jobz='N'; if(eigenvector) jobz='V';
899
900 int_4 n = a.Size(rowa);
901 int_4 lda = a.Step(cola);
902 int_4 info = 0;
903 int_4 lwork = -1;
904 T * work = NULL;
905 T wkget[2];
906
907 b.ReSize(n); b = 0.;
908
909 if (typeid(T) == typeid(r_4) ) {
910 r_4* w = new r_4[n];
911 ssyev_(&jobz,&uplo,&n,(r_4 *)a.Data(),&lda,(r_4 *)w,(r_4 *)wkget,&lwork,&info);
912 lwork = type2i4(&wkget[0],4); /* 3*n-1;*/ work = new T[lwork +GARDMEM];
913 ssyev_(&jobz,&uplo,&n,(r_4 *)a.Data(),&lda,(r_4 *)w,(r_4 *)work,&lwork,&info);
914 if(info==0) for(int i=0;i<n;i++) b(i) = w[i];
915 delete [] w;
916 } else if (typeid(T) == typeid(r_8) ) {
917 r_8* w = new r_8[n];
918 dsyev_(&jobz,&uplo,&n,(r_8 *)a.Data(),&lda,(r_8 *)w,(r_8 *)wkget,&lwork,&info);
919 lwork = type2i4(&wkget[0],8); /* 3*n-1;*/ work = new T[lwork +GARDMEM];
920 dsyev_(&jobz,&uplo,&n,(r_8 *)a.Data(),&lda,(r_8 *)w,(r_8 *)work,&lwork,&info);
921 if(info==0) for(int i=0;i<n;i++) b(i) = w[i];
922 delete [] w;
923 } else if (typeid(T) == typeid(complex<r_4>) ) {
924 r_4* rwork = new r_4[3*n-2 +GARDMEM]; r_4* w = new r_4[n];
925 cheev_(&jobz,&uplo,&n,(complex<r_4> *)a.Data(),&lda,(r_4 *)w
926 ,(complex<r_4> *)wkget,&lwork,(r_4 *)rwork,&info);
927 lwork = type2i4(&wkget[0],4); /* 2*n-1;*/ work = new T[lwork +GARDMEM];
928 cheev_(&jobz,&uplo,&n,(complex<r_4> *)a.Data(),&lda,(r_4 *)w
929 ,(complex<r_4> *)work,&lwork,(r_4 *)rwork,&info);
930 if(info==0) for(int i=0;i<n;i++) b(i) = w[i];
931 delete [] rwork; delete [] w;
932 } else if (typeid(T) == typeid(complex<r_8>) ) {
933 r_8* rwork = new r_8[3*n-2 +GARDMEM]; r_8* w = new r_8[n];
934 zheev_(&jobz,&uplo,&n,(complex<r_8> *)a.Data(),&lda,(r_8 *)w
935 ,(complex<r_8> *)wkget,&lwork,(r_8 *)rwork,&info);
936 lwork = type2i4(&wkget[0],8); /* 2*n-1;*/ work = new T[lwork +GARDMEM];
937 zheev_(&jobz,&uplo,&n,(complex<r_8> *)a.Data(),&lda,(r_8 *)w
938 ,(complex<r_8> *)work,&lwork,(r_8 *)rwork,&info);
939 if(info==0) for(int i=0;i<n;i++) b(i) = w[i];
940 delete [] rwork; delete [] w;
941 } else {
942 if(work) delete [] work; work=NULL;
943 string tn = typeid(T).name();
944 cerr << " LapackServer::LapackEigenSym(a,b) - Unsupported DataType T = " << tn << endl;
945 throw TypeMismatchExc("LapackServer::LapackEigenSym(a,b) - Unsupported DataType (T)");
946 }
947
948 if(work) delete [] work;
949 return(info);
950}
951
952////////////////////////////////////////////////////////////////////////////////////
953/*! Computes the eigen values and eigen vectors of a general squared matrix \b a.
954 Input arrays should have FortranMemoryMapping (column packed).
955 \param a : input general n-by-n matrix
956 \param eval : Vector of eigenvalues (complex double precision)
957 \param evec : Matrix of eigenvector (same order than eval, one vector = one column)
958 \param eigenvector : if true compute (right) eigenvectors, if not only eigenvalues
959 \param a : on return array of eigenvectors
960 \return : return code from lapack driver
961 \verbatim
962 eval : contains the computed eigenvalues.
963 --- For real matrices "a" :
964 Complex conjugate pairs of eigenvalues appear consecutively
965 with the eigenvalue having the positive imaginary part first.
966 evec : the right eigenvectors v(j) are stored one after another
967 in the columns of evec, in the same order as their eigenvalues.
968 --- For real matrices "a" :
969 If the j-th eigenvalue is real, then v(j) = evec(:,j),
970 the j-th column of evec.
971 If the j-th and (j+1)-st eigenvalues form a complex
972 conjugate pair, then v(j) = evec(:,j) + i*evec(:,j+1) and
973 v(j+1) = evec(:,j) - i*evec(:,j+1).
974 \endverbatim
975*/
976
977template <class T>
978int LapackServer<T>::LapackEigen(TArray<T>& a, TVector< complex<r_8> >& eval, TMatrix<T>& evec, bool eigenvector)
979{
980 if ( a.NbDimensions() != 2 )
981 throw(SzMismatchError("LapackServer::LapackEigen(a,b) a NbDimensions() != 2"));
982 int_4 rowa = a.RowsKA();
983 int_4 cola = a.ColsKA();
984 if ( a.Size(rowa) != a.Size(cola))
985 throw(SzMismatchError("LapackServer::LapackEigen(a,b) a Not a square Array"));
986 if (!a.IsPacked(rowa))
987 throw(SzMismatchError("LapackServer::LapackEigen(a,b) a Not Column Packed"));
988
989 char jobvl = 'N';
990 char jobvr = 'N'; if(eigenvector) jobvr='V';
991
992 int_4 n = a.Size(rowa);
993 int_4 lda = a.Step(cola);
994 int_4 info = 0;
995
996 eval.ReSize(n); eval = complex<r_8>(0.,0.);
997 if(eigenvector) {evec.ReSize(n,n); evec = (T) 0.;}
998 int_4 ldvr = n, ldvl = 1;
999
1000 int_4 lwork = -1;
1001 T * work = NULL;
1002 T wkget[2];
1003
1004 if (typeid(T) == typeid(r_4) ) {
1005 r_4* wr = new r_4[n]; r_4* wi = new r_4[n]; r_4* vl = NULL;
1006 sgeev_(&jobvl,&jobvr,&n,(r_4 *)a.Data(),&lda,(r_4 *)wr,(r_4 *)wi,
1007 (r_4 *)vl,&ldvl,(r_4 *)evec.Data(),&ldvr,
1008 (r_4 *)wkget,&lwork,&info);
1009 lwork = type2i4(&wkget[0],4); /* 4*n;*/ work = new T[lwork +GARDMEM];
1010 sgeev_(&jobvl,&jobvr,&n,(r_4 *)a.Data(),&lda,(r_4 *)wr,(r_4 *)wi,
1011 (r_4 *)vl,&ldvl,(r_4 *)evec.Data(),&ldvr,
1012 (r_4 *)work,&lwork,&info);
1013 if(info==0) for(int i=0;i<n;i++) eval(i) = complex<r_8>(wr[i],wi[i]);
1014 delete [] wr; delete [] wi;
1015 } else if (typeid(T) == typeid(r_8) ) {
1016 r_8* wr = new r_8[n]; r_8* wi = new r_8[n]; r_8* vl = NULL;
1017 dgeev_(&jobvl,&jobvr,&n,(r_8 *)a.Data(),&lda,(r_8 *)wr,(r_8 *)wi,
1018 (r_8 *)vl,&ldvl,(r_8 *)evec.Data(),&ldvr,
1019 (r_8 *)wkget,&lwork,&info);
1020 lwork = type2i4(&wkget[0],8); /* 4*n;*/ work = new T[lwork +GARDMEM];
1021 dgeev_(&jobvl,&jobvr,&n,(r_8 *)a.Data(),&lda,(r_8 *)wr,(r_8 *)wi,
1022 (r_8 *)vl,&ldvl,(r_8 *)evec.Data(),&ldvr,
1023 (r_8 *)work,&lwork,&info);
1024 if(info==0) for(int i=0;i<n;i++) eval(i) = complex<r_8>(wr[i],wi[i]);
1025 delete [] wr; delete [] wi;
1026 } else if (typeid(T) == typeid(complex<r_4>) ) {
1027 r_4* rwork = new r_4[2*n +GARDMEM]; r_4* vl = NULL; TVector< complex<r_4> > w(n);
1028 cgeev_(&jobvl,&jobvr,&n,(complex<r_4> *)a.Data(),&lda,(complex<r_4> *)w.Data(),
1029 (complex<r_4> *)vl,&ldvl,(complex<r_4> *)evec.Data(),&ldvr,
1030 (complex<r_4> *)wkget,&lwork,(r_4 *)rwork,&info);
1031 lwork = type2i4(&wkget[0],4); /* 2*n;*/ work = new T[lwork +GARDMEM];
1032 cgeev_(&jobvl,&jobvr,&n,(complex<r_4> *)a.Data(),&lda,(complex<r_4> *)w.Data(),
1033 (complex<r_4> *)vl,&ldvl,(complex<r_4> *)evec.Data(),&ldvr,
1034 (complex<r_4> *)work,&lwork,(r_4 *)rwork,&info);
1035 if(info==0) for(int i=0;i<n;i++) eval(i) = w(i);
1036 delete [] rwork;
1037 } else if (typeid(T) == typeid(complex<r_8>) ) {
1038 r_8* rwork = new r_8[2*n +GARDMEM]; r_8* vl = NULL;
1039 zgeev_(&jobvl,&jobvr,&n,(complex<r_8> *)a.Data(),&lda,(complex<r_8> *)eval.Data(),
1040 (complex<r_8> *)vl,&ldvl,(complex<r_8> *)evec.Data(),&ldvr,
1041 (complex<r_8> *)wkget,&lwork,(r_8 *)rwork,&info);
1042 lwork = type2i4(&wkget[0],8); /* 2*n;*/ work = new T[lwork +GARDMEM];
1043 zgeev_(&jobvl,&jobvr,&n,(complex<r_8> *)a.Data(),&lda,(complex<r_8> *)eval.Data(),
1044 (complex<r_8> *)vl,&ldvl,(complex<r_8> *)evec.Data(),&ldvr,
1045 (complex<r_8> *)work,&lwork,(r_8 *)rwork,&info);
1046 delete [] rwork;
1047 } else {
1048 if(work) delete [] work; work=NULL;
1049 string tn = typeid(T).name();
1050 cerr << " LapackServer::LapackEigen(a,b) - Unsupported DataType T = " << tn << endl;
1051 throw TypeMismatchExc("LapackServer::LapackEigen(a,b) - Unsupported DataType (T)");
1052 }
1053
1054 if(work) delete [] work;
1055 return(info);
1056}
1057
1058
1059
1060
1061///////////////////////////////////////////////////////////////
1062#ifdef __CXX_PRAGMA_TEMPLATES__
1063#pragma define_template LapackServer<r_4>
1064#pragma define_template LapackServer<r_8>
1065#pragma define_template LapackServer< complex<r_4> >
1066#pragma define_template LapackServer< complex<r_8> >
1067#endif
1068
1069#if defined(ANSI_TEMPLATES) || defined(GNU_TEMPLATES)
1070template class LapackServer<r_4>;
1071template class LapackServer<r_8>;
1072template class LapackServer< complex<r_4> >;
1073template class LapackServer< complex<r_8> >;
1074#endif
1075
1076#if defined(OS_LINUX)
1077// Pour le link avec f2c sous Linux
1078extern "C" {
1079 void MAIN__();
1080}
1081
1082void MAIN__()
1083{
1084 cerr << "MAIN__() function for linking with libf2c.a " << endl;
1085 cerr << " This function should never be called !!! " << endl;
1086 throw PError("MAIN__() should not be called - see intflapack.cc");
1087}
1088#endif
Note: See TracBrowser for help on using the repository browser.