1 //
2 // LAPACK++ 1.1 Linear Algebra Package 1.1
3 // University of Tennessee, Knoxvilee, TN.
4 // Oak Ridge National Laboratory, Oak Ridge, TN.
5 // Authors: J. J. Dongarra, E. Greaser, R. Pozo, D. Walker
6 // (C) 1992-1996 All Rights Reserved
7 //
8 // NOTICE
9 //
10 // Permission to use, copy, modify, and distribute this software and
11 // its documentation for any purpose and without fee is hereby granted
12 // provided that the above copyright notice appear in all copies and
13 // that both the copyright notice and this permission notice appear in
14 // supporting documentation.
15 //
16 // Neither the Institutions (University of Tennessee, and Oak Ridge National
17 // Laboratory) nor the Authors make any representations about the suitability
18 // of this software for any purpose. This software is provided ``as is''
19 // without express or implied warranty.
20 //
21 // LAPACK++ was funded in part by the U.S. Department of Energy, the
22 // National Science Foundation and the State of Tennessee.
23
24 #ifdef HAVE_CONFIG_H
25 # include <config.h>
26 #endif
27
28 #include <iostream>
29 #include "lapack.h"
30 #include "lapackc.h"
31 #include "lafnames.h"
32 #include LA_EXCEPTION_H
33 #include LA_GEN_MAT_DOUBLE_H
34 #include LA_VECTOR_DOUBLE_H
35 #include LA_VECTOR_LONG_INT_H
36 #ifdef LA_COMPLEX_SUPPORT
37 # include LA_GEN_MAT_COMPLEX_H
38 # include LA_VECTOR_COMPLEX_H
39 #endif
40 #include LA_SPD_MAT_DOUBLE_H
41 #include LA_SYMM_MAT_DOUBLE_H
42 #include "blas3pp.h"
43
44 #include LA_SOLVE_DOUBLE_H
45 #include LA_UTIL_H
46
47
LaLinearSolve(const LaGenMatDouble & A,LaGenMatDouble & X,const LaGenMatDouble & B)48 void LaLinearSolve( const LaGenMatDouble& A, LaGenMatDouble& X,
49 const LaGenMatDouble& B)
50 {
51 int M = A.size(0), N = A.size(1);
52
53 if ( M == N )
54 LaLULinearSolve(A, X, B);
55 else
56 LaQRLinearSolve(A, X, B);
57 }
58
LaLinearSolve(const LaSpdMatDouble & A,LaGenMatDouble & X,LaGenMatDouble & B)59 void LaLinearSolve(const LaSpdMatDouble &A, LaGenMatDouble& X,
60 LaGenMatDouble& B )
61 {
62 LaCholLinearSolve(A, X, B );
63 }
64
LaLinearSolve(const LaSymmMatDouble & A,LaGenMatDouble & X,const LaGenMatDouble & B)65 void LaLinearSolve(const LaSymmMatDouble &A, LaGenMatDouble& X,
66 const LaGenMatDouble& B )
67 {
68 LaCholLinearSolve(A, X, B );
69 }
70
LaLinearSolveIP(LaSpdMatDouble & A,LaGenMatDouble & X,LaGenMatDouble & B)71 void LaLinearSolveIP(LaSpdMatDouble &A, LaGenMatDouble& X, LaGenMatDouble& B )
72 {
73 LaCholLinearSolveIP(A, X, B );
74 }
75
LaLinearSolveIP(LaSymmMatDouble & A,LaGenMatDouble & X,const LaGenMatDouble & B)76 void LaLinearSolveIP(LaSymmMatDouble &A, LaGenMatDouble& X,
77 const LaGenMatDouble& B )
78 {
79 LaCholLinearSolveIP(A, X, B );
80 }
81
LaLinearSolveIP(LaGenMatDouble & A,LaGenMatDouble & X,const LaGenMatDouble & B)82 void LaLinearSolveIP( LaGenMatDouble& A, LaGenMatDouble& X,
83 const LaGenMatDouble& B)
84 {
85 int M = A.size(0), N = A.size(1);
86
87 if ( M == N )
88 LaLULinearSolveIP(A, X, B);
89 else
90 LaQRLinearSolveIP(A, X, B);
91 }
92
LaLULinearSolve(const LaGenMatDouble & A,LaGenMatDouble & X,const LaGenMatDouble & B)93 void LaLULinearSolve(const LaGenMatDouble& A, LaGenMatDouble& X,
94 const LaGenMatDouble& B )
95 {
96 LaGenMatDouble A1(A); // exception if out of memory
97 LaLULinearSolveIP(A1, X, B);
98 }
99
100
101 // General LU Solver
102 //
103 // N x N N x nrhs N x nrhs
104 //
LaLULinearSolveIP(LaGenMatDouble & A,LaGenMatDouble & X,const LaGenMatDouble & B)105 void LaLULinearSolveIP( LaGenMatDouble& A, LaGenMatDouble& X,
106 const LaGenMatDouble& B )
107 {
108 #ifndef HPPA
109 const char fname[] = "LaLULinearSolveIP(LaGenMatDouble &A, &X, &B)";
110 #else
111 char *fname = NULL; // HP C++ does not support string initalization!
112 #endif
113
114
115 // let's not worry about non-unit column strides for the moment
116 if ( A.inc(0) != 1 || A.inc(1) != 1)
117 throw(LaException(fname, "A is non-contiguous."));
118
119 if (!(X.size(0) == B.size(0) && X.size(1) == B.size(1)))
120 throw(LaException(fname, "X and B are non-conformant."));
121
122 X.inject(B); // will throw exception if not conformant
123
124
125 // in the future this can call the linear least square routines
126 // to handle non-square matrices
127
128 if (A.size(0) != A.size(1))
129 throw(LaException(fname, "Square matrix expected.\n"));
130
131 if (A.size(1) != X.size(0))
132 throw(LaException(fname, "A and X are non-comformant."));
133
134 integer info = 0;
135 int M = A.size(0);
136 integer Ml = M;
137 //integer N = A.size(1);
138 integer K = X.size(1);
139 integer lda = A.inc(0) * A.gdim(0);
140 integer ldx = X.inc(0) * X.gdim(0);
141
142 LaVectorLongInt ipiv( M);
143
144 F77NAME(dgesv) (&Ml, &K, &A(0, 0), &lda, &ipiv(0), &X(0, 0), &ldx, &info);
145
146 if (info < 0)
147 throw(LaException(fname, "Internal error in LAPACK: DGESV() with illegal argument value"));
148 else if (info > 0)
149 throw(LaException(fname, "Internal error in LAPACK: DGESV() Factor U was exactly singular"));
150 }
151
152
153
154
155
LaQRLinearSolve(const LaGenMatDouble & A,LaGenMatDouble & X,const LaGenMatDouble & B)156 void LaQRLinearSolve(const LaGenMatDouble& A, LaGenMatDouble& X,
157 const LaGenMatDouble& B )
158 {
159 LaGenMatDouble A1(A);
160 LaQRLinearSolveIP(A1, X, B);
161 }
162
163
164 // General QR solver
165 //
166 // M x N N x nrhs M x nrhs
167 //
LaQRLinearSolveIP(LaGenMatDouble & A,LaGenMatDouble & X,const LaGenMatDouble & B)168 void LaQRLinearSolveIP(LaGenMatDouble& A, LaGenMatDouble& X,
169 const LaGenMatDouble& B )
170 {
171 #ifndef HPPA
172 const char fname[] = "LaQRLinearSolveIP(LaGenMatDouble &A, &X, &B)";
173 #else
174 char *fname = NULL; // HP C++ does not support string initalization!
175 #endif
176
177 // let's not worry about non-unit column strides for the moment
178 if ( A.inc(0) != 1 || A.inc(1) != 1)
179 throw(LaException(fname, "A is non-contiguous."));
180 if ( A.size(0) == 0 || A.size(1) == 0 )
181 throw(LaException(fname, "Matrix A is empty; one dimension is zero."));
182
183 if (!( A.size(0) == B.size(0) &&
184 A.size(1) == X.size(0) &&
185 X.size(1) == B.size(1) ))
186 throw(LaException(fname, "input matrices are non-conformant."));
187
188 integer info = 0;
189 int M = A.size(0);
190 int N = A.size(1);
191 integer Ml = M;
192 integer Nl = N;
193 int nrhs = X.size(1);
194 integer nrhsl = nrhs;
195 integer lda = A.inc(0) * A.gdim(0);
196
197 int nb = LaEnvBlockSize("DGELS", A);
198 integer lwork = M * N + nb * std::max(M * N, nrhs);
199 //std::cout << fname << ": nb= " << nb << " lwork=" << lwork << std::endl;
200
201 LaVectorDouble WORK(lwork);
202 char trans = 'N';
203
204 if (M != N)
205 {
206 // Typically is A non-square, so we need to create tmp X because
207 // X is N x nrhs, while B is M x nrhs. We need to make copies of
208 // these so that the routine won't corrupt data around X and B.
209
210 LaGenMatDouble Xtmp(std::max(M, N), nrhs);
211 integer ldx = Xtmp.inc(0) * Xtmp.gdim(0);
212
213 // Copy B into the temporary X matrix that is passed to dgels()
214 Xtmp(LaIndex(0, M - 1), LaIndex()).inject( B );
215
216 F77NAME(dgels) (&trans, &Ml, &Nl, &nrhsl, &A(0, 0), &lda, &Xtmp(0, 0),
217 &ldx, &WORK(0), &lwork, &info);
218
219 // And copy the result from the larger matrix back into
220 // the actual result matrix.
221 X.inject(Xtmp(LaIndex(0, N - 1), LaIndex()));
222 }
223 else
224 {
225 integer ldx = X.inc(0) * X.gdim(0);
226
227 // Copy B into the X matrix that is passed to dgels()
228 X.inject( B );
229
230 F77NAME(dgels) (&trans, &Ml, &Nl, &nrhsl, &A(0, 0), &lda, &X(0, 0),
231 &ldx, &WORK(0), &lwork, &info);
232 }
233
234
235 // this shouldn't really happen.
236 //
237 if (info < 0)
238 throw(LaException(fname, "Internal error in LAPACK: SGELS()"));
239
240 }
241 // ////////////////////////////////////////////////////////////
242
243 #ifdef LA_COMPLEX_SUPPORT
LaLinearSolve(const LaGenMatComplex & A,LaGenMatComplex & X,const LaGenMatComplex & B)244 void LaLinearSolve( const LaGenMatComplex& A, LaGenMatComplex& X,
245 const LaGenMatComplex& B)
246 {
247 int M = A.size(0), N = A.size(1);
248
249 if ( M == N )
250 LaLULinearSolve(A, X, B);
251 else
252 LaQRLinearSolve(A, X, B);
253 }
254
LaLinearSolveIP(LaGenMatComplex & A,LaGenMatComplex & X,const LaGenMatComplex & B)255 void LaLinearSolveIP( LaGenMatComplex& A, LaGenMatComplex& X,
256 const LaGenMatComplex& B)
257 {
258 int M = A.size(0), N = A.size(1);
259
260 if ( M == N )
261 LaLULinearSolveIP(A, X, B);
262 else
263 LaQRLinearSolveIP(A, X, B);
264 }
265
LaLULinearSolve(const LaGenMatComplex & A,LaGenMatComplex & X,const LaGenMatComplex & B)266 void LaLULinearSolve(const LaGenMatComplex& A, LaGenMatComplex& X,
267 const LaGenMatComplex& B )
268 {
269 LaGenMatComplex A1(A); // exception if out of memory
270 LaLULinearSolveIP(A1, X, B);
271 }
272
273
274 // General LU Solver
275 //
276 // N x N N x nrhs N x nrhs
277 //
LaLULinearSolveIP(LaGenMatComplex & A,LaGenMatComplex & X,const LaGenMatComplex & B)278 void LaLULinearSolveIP( LaGenMatComplex& A, LaGenMatComplex& X,
279 const LaGenMatComplex& B )
280 {
281 #ifndef HPPA
282 const char fname[] = "LaLULinearSolveIP(LaGenMatComplex &A, &X, &B)";
283 #else
284 char *fname = NULL; // HP C++ does not support string initalization!
285 #endif
286
287
288 // let's not worry about non-unit column strides for the moment
289 if ( A.inc(0) != 1 || A.inc(1) != 1)
290 throw(LaException(fname, "A is non-contiguous."));
291
292 if (!(X.size(0) == B.size(0) && X.size(1) == B.size(1)))
293 throw(LaException(fname, "X and B are non-conformant."));
294
295 X.inject(B); // will throw exception if not conformant
296
297
298 // in the future this can call the linear least square routines
299 // to handle non-square matrices
300
301 if (A.size(0) != A.size(1))
302 throw(LaException(fname, "Square matrix expected.\n"));
303
304 if (A.size(1) != X.size(0))
305 throw(LaException(fname, "A and X are non-comformant."));
306
307 integer info = 0;
308 int M = A.size(0);
309 integer Ml = M;
310 //integer N = A.size(1);
311 integer K = X.size(1);
312 integer lda = A.inc(0) * A.gdim(0);
313 integer ldx = X.inc(0) * X.gdim(0);
314
315 LaVectorLongInt ipiv( M);
316
317 F77NAME(zgesv) (&Ml, &K, &A(0, 0), &lda, &ipiv(0), &X(0, 0), &ldx, &info);
318
319 if (info < 0)
320 throw(LaException(fname, "Internal error in LAPACK: DGESV() with illegal argument value"));
321 else if (info > 0)
322 throw(LaException(fname, "Internal error in LAPACK: DGESV() Factor U was exactly singular"));
323 }
324
325
326
327
328
LaQRLinearSolve(const LaGenMatComplex & A,LaGenMatComplex & X,const LaGenMatComplex & B)329 void LaQRLinearSolve(const LaGenMatComplex& A, LaGenMatComplex& X,
330 const LaGenMatComplex& B )
331 {
332 LaGenMatComplex A1(A);
333 LaQRLinearSolveIP(A1, X, B);
334 }
335
336
337 // General QR solver
338 //
339 // M x N N x nrhs M x nrhs
340 //
LaQRLinearSolveIP(LaGenMatComplex & A,LaGenMatComplex & X,const LaGenMatComplex & B)341 void LaQRLinearSolveIP(LaGenMatComplex& A, LaGenMatComplex& X, const LaGenMatComplex& B )
342 {
343 #ifndef HPPA
344 const char fname[] = "LaQRLinearSolveIP(LaGenMatComplex &A, &X, &B)";
345 #else
346 char *fname = NULL; // HP C++ does not support string initalization!
347 #endif
348
349 // let's not worry about non-unit column strides for the moment
350 if ( A.inc(0) != 1 || A.inc(1) != 1)
351 throw(LaException(fname, "A is non-contiguous."));
352 if ( A.size(0) == 0 || A.size(1) == 0 )
353 throw(LaException(fname, "Matrix A is empty; one dimension is zero."));
354
355 if (!( A.size(0) == B.size(0) &&
356 A.size(1) == X.size(0) &&
357 X.size(1) == B.size(1) ))
358 throw(LaException(fname, "input matrices are non-conformant."));
359
360 integer info = 0;
361 int M = A.size(0);
362 int N = A.size(1);
363 integer Ml = M;
364 integer Nl = N;
365 int nrhs = X.size(1);
366 integer nrhsl = nrhs;
367 integer lda = A.inc(0) * A.gdim(0);
368
369
370 //int nb = 32;
371 int nb = LaEnvBlockSize("ZGELS", A);
372 integer lwork = M * N + nb * std::max(M * N, nrhs);
373 //std::cout << "Block size: " << nb << std::endl;
374
375 LaVectorComplex WORK(lwork);
376 char trans = 'N';
377
378 if (M != N)
379 {
380 // Typically is A non-square, so we need to create tmp X because
381 // X is N x nrhs, while B is M x nrhs. We need to make copies of
382 // these so that the routine won't corrupt data around X and B.
383
384 LaGenMatComplex Xtmp(std::max(M, N), nrhs);
385 integer ldx = Xtmp.inc(0) * Xtmp.gdim(0);
386
387 // Copy B into the temporary X matrix that is passed to zgels()
388 Xtmp(LaIndex(0, M - 1), LaIndex()).inject( B );
389
390 F77NAME(zgels) (&trans, &Ml, &Nl, &nrhsl, &A(0, 0), &lda, &Xtmp(0, 0),
391 &ldx, &WORK(0), &lwork, &info);
392
393 // And copy the result from the larger matrix back into
394 // the actual result matrix.
395 X.inject(Xtmp(LaIndex(0, N - 1), LaIndex()));
396 }
397 else
398 {
399 integer ldx = X.inc(0) * X.gdim(0);
400
401 // Copy B into the X matrix that is passed to dgels()
402 X.inject( B );
403
404 F77NAME(zgels) (&trans, &Ml, &Nl, &nrhsl, &A(0, 0), &lda, &X(0, 0),
405 &ldx, &WORK(0), &lwork, &info);
406 }
407
408
409 // this shouldn't really happen.
410 //
411 if (info < 0)
412 throw(LaException(fname, "Internal error in LAPACK: ZGELS()"));
413
414 }
415
416 #endif // LA_COMPLEX_SUPPORT
417
418
419
420 // ////////////////////////////////////////////////////////////
421
LaCholLinearSolve(const LaSpdMatDouble & A,LaGenMatDouble & X,LaGenMatDouble & B)422 void LaCholLinearSolve( const LaSpdMatDouble& A, LaGenMatDouble& X,
423 LaGenMatDouble& B )
424 {
425 LaSpdMatDouble A1(A);
426 LaCholLinearSolveIP(A1, X, B);
427 }
428
LaCholLinearSolve(const LaSymmMatDouble & A,LaGenMatDouble & X,const LaGenMatDouble & B)429 void LaCholLinearSolve( const LaSymmMatDouble& A, LaGenMatDouble& X,
430 const LaGenMatDouble& B )
431 {
432 LaSymmMatDouble A1(A);
433 LaCholLinearSolveIP(A1, X, B);
434 }
435
436 // A is NxN, X is MxN and B is MxN
437 //
LaCholLinearSolveIP(LaSpdMatDouble & A,LaGenMatDouble & X,LaGenMatDouble & B)438 void LaCholLinearSolveIP( LaSpdMatDouble& A, LaGenMatDouble& X,
439 LaGenMatDouble& B )
440 {
441 #ifndef HPPA
442 const char fname[] = "LaCholLinearSolveIP(LaSpdMatDouble &A, &X, &B)";
443 #else
444 char *fname = NULL; // HP C++ does not support string initalization!
445 #endif
446
447 // let's not worry about non-unit column strides for the moment
448 if ( A.inc(0) != 1 || A.inc(1) != 1)
449 throw(LaException(fname, "A is non-contiguous."));
450
451 if (!(X.size(0) == B.size(0) && X.size(1) == B.size(1)))
452 throw(LaException(fname, "X and B are non-conformant."));
453
454 X.inject(B); // will throw exception if not conformant
455
456
457 // in the future this can call the linear least square routines
458 // to handle non-square matrices
459
460 if (A.size(0) != A.size(1))
461 throw(LaException(fname, "Square matrix expected.\n"));
462
463 if (A.size(1) != X.size(0))
464 throw(LaException(fname, "A and X are non-comformant."));
465
466 integer info = 0;
467 integer M = A.size(0);
468 //integer N = A.size(1);
469 integer K = X.size(1);
470 integer lda = A.inc(0) * A.gdim(0);
471 integer ldx = X.inc(0) * X.gdim(0);
472 char uplo = 'L';
473
474 F77NAME(dposv) (&uplo, &M, &K, &A(0, 0), &lda, &X(0, 0), &ldx, &info);
475
476 // this shouldn't really happen.
477 //
478 if (info < 0)
479 throw(LaException(fname, "Internal error in LAPACK: SGESV()"));
480
481 if (info > 0)
482 throw (LaException(fname, "A is not symmetric-positive-definite."));
483
484
485 }
486
LaCholLinearSolveIP(LaSymmMatDouble & A,LaGenMatDouble & X,const LaGenMatDouble & B)487 void LaCholLinearSolveIP( LaSymmMatDouble& A, LaGenMatDouble& X,
488 const LaGenMatDouble& B )
489 {
490 #ifndef HPPA
491 const char fname[] = "LaCholLinearSolveIP(LaSymmMatDouble &A, &X, &B)";
492 #else
493 char *fname = NULL; // HP C++ does not support string initalization!
494 #endif
495
496 // let's not worry about non-unit column strides for the moment
497 if ( A.inc(0) != 1 || A.inc(1) != 1)
498 throw(LaException(fname, "A is non-contiguous."));
499
500 if (!(X.size(0) == B.size(0) && X.size(1) == B.size(1)))
501 throw(LaException(fname, "X and B are non-conformant."));
502
503 X.inject(B); // will throw exception if not conformant
504
505
506 // in the future this can call the linear least square routines
507 // to handle non-square matrices
508
509 if (A.size(0) != A.size(1))
510 throw(LaException(fname, "Square matrix expected.\n"));
511
512 if (A.size(1) != X.size(0))
513 throw(LaException(fname, "A and X are non-comformant."));
514
515 integer info = 0;
516 integer M = A.size(0);
517 //integer N = A.size(1);
518 integer K = X.size(1);
519 integer lda = A.inc(0) * A.gdim(0);
520 integer ldx = X.inc(0) * X.gdim(0);
521 char uplo = 'L';
522
523 LaVectorLongInt ipiv(M);
524 integer lwork = -1;
525 LaVectorDouble work(1);
526 // Workspace query
527 F77NAME(dsysv) (&uplo, &M, &K, &A(0, 0), &lda, &ipiv(0), &X(0, 0), &ldx,
528 &work(0), &lwork, &info);
529 lwork = integer(work(0));
530 work.resize(lwork, 1);
531
532 F77NAME(dsysv) (&uplo, &M, &K, &A(0, 0), &lda, &ipiv(0), &X(0, 0), &ldx,
533 &work(0), &lwork, &info);
534
535 // this shouldn't really happen.
536 //
537 if (info < 0)
538 throw(LaException(fname, "Internal error in LAPACK: DSYSV()"));
539
540 if (info > 0)
541 throw (LaException(fname, "Matrix is singular."));
542
543
544 }
545
546
LUFactorizeIP(LaGenMatDouble & GM,LaVectorLongInt & PIV)547 void LUFactorizeIP(LaGenMatDouble &GM, LaVectorLongInt &PIV)
548 {
549 integer m = GM.size(0), n = GM.size(1), lda = GM.gdim(0);
550 integer info = 0;
551 assert(PIV.size() >= (m < n ? m : n));
552
553 // Copied from LaGenMatFactorize in fmd.h of LAPACK++, but the
554 // version there didn't admit to being in-place modifying.
555 F77NAME(dgetrf)(&m, &n, &GM(0, 0), &lda, &(PIV(0)), &info);
556
557 if (info < 0)
558 throw LaException("LUFactorizeIP", "Error in argument");
559 }
560
561 #ifdef LA_COMPLEX_SUPPORT
LUFactorizeIP(LaGenMatComplex & GM,LaVectorLongInt & PIV)562 void LUFactorizeIP(LaGenMatComplex &GM, LaVectorLongInt &PIV)
563 {
564 integer m = GM.size(0), n = GM.size(1), lda = GM.gdim(0);
565 integer info = 0;
566 assert(PIV.size() >= (m < n ? m : n));
567
568 // Copied from LaGenMatFactorize in fmd.h of LAPACK++, but the
569 // version there didn't admit to being in-place modifying.
570 F77NAME(zgetrf)(&m, &n, &GM(0, 0), &lda, &(PIV(0)), &info);
571
572 if (info < 0)
573 throw LaException("LUFactorizeIP", "Error in argument");
574 }
575
LaLUInverseIP(LaGenMatComplex & A,LaVectorLongInt & PIV)576 void LaLUInverseIP(LaGenMatComplex &A, LaVectorLongInt &PIV)
577 {
578 LaVectorComplex work; // will be resized in other function
579 LaLUInverseIP(A, PIV, work);
580 }
LaLUInverseIP(LaGenMatComplex & A,LaVectorLongInt & PIV,LaVectorComplex & work)581 void LaLUInverseIP(LaGenMatComplex &A, LaVectorLongInt &PIV, LaVectorComplex &work)
582 {
583 integer N = A.size(1), lda = A.gdim(0), info = 0;
584 if (A.size(0) != A.size(1))
585 throw LaException("LaLUInverseIP", "Input must be square");
586 integer W = work.size();
587
588 // Check for minimum work size
589 if ( W < A.size(0) )
590 {
591 int nb = LaEnvBlockSize("ZGETRI", A);
592 W = N * nb;
593 work.resize(W, 1);
594 }
595
596 F77NAME(zgetri)(&N, &(A(0, 0)), &lda, &(PIV(0)), &work(0), &W, &info);
597 if (info < 0)
598 throw LaException("LaLUInverseIP", "Error in zgetri argument");
599 if (info > 0)
600 throw LaException("LaLuInverseIP", "Matrix is singlular, cannot compute inverse");
601 }
602 #endif // LA_COMPLEX_SUPPORT
603
LaLUInverseIP(LaGenMatDouble & A,LaVectorLongInt & PIV)604 void LaLUInverseIP(LaGenMatDouble &A, LaVectorLongInt &PIV)
605 {
606 LaVectorDouble work; // will be resized in other function
607 LaLUInverseIP(A, PIV, work);
608 #if 0
609 // above code is shorter - remove this code soon.
610 integer N = A.size(1), lda = A.gdim(0), info = 0;
611 if (A.size(0) != A.size(1))
612 throw LaException("LaLUInverseIP", "Input must be square");
613 int nb = LaEnvBlockSize("DGETRI", A);
614 integer W = N * nb;
615 LaVectorDouble work(W);
616
617 F77NAME(dgetri)(&N, &(A(0, 0)), &lda, &(PIV(0)), &work(0), &W, &info);
618 if (info < 0)
619 throw LaException("inv", "Error in dgetri argument");
620 #endif
621 }
622
LaLUInverseIP(LaGenMatDouble & A,LaVectorLongInt & PIV,LaVectorDouble & work)623 void LaLUInverseIP(LaGenMatDouble &A, LaVectorLongInt &PIV, LaVectorDouble &work)
624 {
625 integer N = A.size(1), lda = A.gdim(0), info = 0;
626 if (A.size(0) != A.size(1))
627 throw LaException("LaLUInverseIP", "Input must be square");
628 integer W = work.size();
629
630 // Check for minimum work size
631 if ( W < A.size(0) )
632 {
633 int nb = LaEnvBlockSize("DGETRI", A);
634 W = N * nb;
635 work.resize(W, 1);
636 }
637
638 F77NAME(dgetri)(&N, &(A(0, 0)), &lda, &(PIV(0)), &work(0), &W, &info);
639 if (info < 0)
640 throw LaException("LaLUInverseIP", "Error in dgetri argument");
641 if (info > 0)
642 throw LaException("LaLuInverseIP", "Matrix is singlular, cannot compute inverse");
643 }
644