001 // --- BEGIN LICENSE BLOCK ---
002 /*
003 * Copyright (c) 2009, Mikio L. Braun
004 * All rights reserved.
005 *
006 * Redistribution and use in source and binary forms, with or without
007 * modification, are permitted provided that the following conditions are
008 * met:
009 *
010 * * Redistributions of source code must retain the above copyright
011 * notice, this list of conditions and the following disclaimer.
012 *
013 * * Redistributions in binary form must reproduce the above
014 * copyright notice, this list of conditions and the following
015 * disclaimer in the documentation and/or other materials provided
016 * with the distribution.
017 *
018 * * Neither the name of the Technische Universit??t Berlin nor the
019 * names of its contributors may be used to endorse or promote
020 * products derived from this software without specific prior
021 * written permission.
022 *
023 * THIS SOFTWARE IS PROVIDED BY THE COPYRIGHT HOLDERS AND CONTRIBUTORS
024 * "AS IS" AND ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT
025 * LIMITED TO, THE IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR
026 * A PARTICULAR PURPOSE ARE DISCLAIMED. IN NO EVENT SHALL THE COPYRIGHT
027 * HOLDER OR CONTRIBUTORS BE LIABLE FOR ANY DIRECT, INDIRECT, INCIDENTAL,
028 * SPECIAL, EXEMPLARY, OR CONSEQUENTIAL DAMAGES (INCLUDING, BUT NOT
029 * LIMITED TO, PROCUREMENT OF SUBSTITUTE GOODS OR SERVICES; LOSS OF USE,
030 * DATA, OR PROFITS; OR BUSINESS INTERRUPTION) HOWEVER CAUSED AND ON ANY
031 * THEORY OF LIABILITY, WHETHER IN CONTRACT, STRICT LIABILITY, OR TORT
032 * (INCLUDING NEGLIGENCE OR OTHERWISE) ARISING IN ANY WAY OUT OF THE USE
033 * OF THIS SOFTWARE, EVEN IF ADVISED OF THE POSSIBILITY OF SUCH DAMAGE.
034 */
035 // --- END LICENSE BLOCK ---
036
037 package org.jblas.la;
038
039 import org.jblas.la.exceptions.LapackException;
040 import org.jblas.la.exceptions.LapackArgumentException;
041 import org.jblas.core.ComplexDouble;
042 import org.jblas.core.ComplexFloat;
043 import org.jblas.la.exceptions.LapackConvergenceException;
044
045 //import edu.ida.core.OutputValue;
046
047 /**
048 * This class provides a cleaner direct interface to the BLAS routines by
049 * extracting the parameters of the matrices from the matrices itself.
050 *
051 * For example, you can just pass the vector and do not have to pass the length,
052 * corresponding DoubleBuffer, offset and step size explicitly.
053 *
054 * Currently, all the general matrix routines are implemented.
055 *
056 */
057 public class SimpleBlas {
058 /***************************************************************************
059 * BLAS Level 1
060 */
061
062 /** Compute x <-> y (swap two matrices) */
063 public static DoubleMatrix swap(DoubleMatrix x, DoubleMatrix y) {
064 //Blas.dswap(x.length, x.data, 0, 1, y.data, 0, 1);
065 JavaBlas.rswap(x.length, x.data, 0, 1, y.data, 0, 1);
066 return y;
067 }
068
069 /** Compute x <- alpha * x (scale a matrix) */
070 public static DoubleMatrix scal(double alpha, DoubleMatrix x) {
071 Blas.dscal(x.length, alpha, x.data, 0, 1);
072 return x;
073 }
074
075 public static ComplexDoubleMatrix scal(ComplexDouble alpha, ComplexDoubleMatrix x) {
076 Blas.zscal(x.length, alpha, x.data, 0, 1);
077 return x;
078 }
079
080 /** Compute y <- x (copy a matrix) */
081 public static DoubleMatrix copy(DoubleMatrix x, DoubleMatrix y) {
082 //Blas.dcopy(x.length, x.data, 0, 1, y.data, 0, 1);
083 JavaBlas.rcopy(x.length, x.data, 0, 1, y.data, 0, 1);
084 return y;
085 }
086
087 public static ComplexDoubleMatrix copy(ComplexDoubleMatrix x, ComplexDoubleMatrix y) {
088 Blas.zcopy(x.length, x.data, 0, 1, y.data, 0, 1);
089 return y;
090 }
091
092 /** Compute y <- alpha * x + y (elementwise addition) */
093 public static DoubleMatrix axpy(double da, DoubleMatrix dx, DoubleMatrix dy) {
094 //Blas.daxpy(dx.length, da, dx.data, 0, 1, dy.data, 0, 1);
095 JavaBlas.raxpy(dx.length, da, dx.data, 0, 1, dy.data, 0, 1);
096
097 return dy;
098 }
099
100 public static ComplexDoubleMatrix axpy(ComplexDouble da, ComplexDoubleMatrix dx, ComplexDoubleMatrix dy) {
101 Blas.zaxpy(dx.length, da, dx.data, 0, 1, dy.data, 0, 1);
102 return dy;
103 }
104
105 /** Compute x^T * y (dot product) */
106 public static double dot(DoubleMatrix x, DoubleMatrix y) {
107 //return Blas.ddot(x.length, x.data, 0, 1, y.data, 0, 1);
108 return JavaBlas.rdot(x.length, x.data, 0, 1, y.data, 0, 1);
109 }
110
111 /** Compute x^T * y (dot product) */
112 public static ComplexDouble dotc(ComplexDoubleMatrix x, ComplexDoubleMatrix y) {
113 return Blas.zdotc(x.length, x.data, 0, 1, y.data, 0, 1);
114 }
115
116 /** Compute x^T * y (dot product) */
117 public static ComplexDouble dotu(ComplexDoubleMatrix x, ComplexDoubleMatrix y) {
118 return Blas.zdotu(x.length, x.data, 0, 1, y.data, 0, 1);
119 }
120
121 /** Compute || x ||_2 (2-norm) */
122 public static double nrm2(DoubleMatrix x) {
123 return Blas.dnrm2(x.length, x.data, 0, 1);
124 }
125
126 public static double nrm2(ComplexDoubleMatrix x) {
127 return Blas.dznrm2(x.length, x.data, 0, 1);
128 }
129
130 /** Compute || x ||_1 (1-norm, sum of absolute values) */
131 public static double asum(DoubleMatrix x) {
132 return Blas.dasum(x.length, x.data, 0, 1);
133 }
134
135 public static double asum(ComplexDoubleMatrix x) {
136 return Blas.dzasum(x.length, x.data, 0, 1);
137 }
138
139 /**
140 * Compute index of element with largest absolute value (index of absolute
141 * value maximum)
142 */
143 public static int iamax(DoubleMatrix x) {
144 return Blas.idamax(x.length, x.data, 0, 1) - 1;
145 }
146
147 public static int iamax(ComplexDoubleMatrix x) {
148 return Blas.izamax(x.length, x.data, 0, 1);
149 }
150
151 /***************************************************************************
152 * BLAS Level 2
153 */
154
155 /**
156 * Compute y <- alpha*op(a)*x + beta * y (general matrix vector
157 * multiplication)
158 */
159 public static DoubleMatrix gemv(double alpha, DoubleMatrix a,
160 DoubleMatrix x, double beta, DoubleMatrix y) {
161 if (false) {
162 Blas.dgemv('N', a.rows, a.columns, alpha, a.data, 0, a.rows, x.data, 0,
163 1, beta, y.data, 0, 1);
164 }
165 else {
166 if (beta == 0.0) {
167 for (int i = 0; i < y.length; i++)
168 y.data[i] = 0.0;
169
170 for (int j = 0; j < a.columns; j++)
171 for (int i = 0; i < a.rows; i++)
172 y.data[i] += a.get(i, j) * x.get(j);
173 }
174 else {
175 for (int j = 0; j < a.columns; j++)
176 for (int i = 0; i < a.rows; i++)
177 y.data[j] = a.get(i, j) * x.get(i) + y.data[j];
178 }
179 }
180 return y;
181 }
182
183 /** Compute A <- alpha * x * y^T + A (general rank-1 update) */
184 public static DoubleMatrix ger(double alpha, DoubleMatrix x,
185 DoubleMatrix y, DoubleMatrix a) {
186 Blas.dger(a.rows, a.columns, alpha, x.data, 0, 1, y.data, 0, 1, a.data,
187 0, a.rows);
188 return a;
189 }
190
191 /** Compute A <- alpha * x * y^T + A (general rank-1 update) */
192 public static ComplexDoubleMatrix geru(ComplexDouble alpha, ComplexDoubleMatrix x,
193 ComplexDoubleMatrix y, ComplexDoubleMatrix a) {
194 Blas.zgeru(a.rows, a.columns, alpha, x.data, 0, 1, y.data, 0, 1, a.data,
195 0, a.rows);
196 return a;
197 }
198
199 /** Compute A <- alpha * x * y^H + A (general rank-1 update) */
200 public static ComplexDoubleMatrix gerc(ComplexDouble alpha, ComplexDoubleMatrix x,
201 ComplexDoubleMatrix y, ComplexDoubleMatrix a) {
202 Blas.zgerc(a.rows, a.columns, alpha, x.data, 0, 1, y.data, 0, 1, a.data,
203 0, a.rows);
204 return a;
205 }
206
207 /***************************************************************************
208 * BLAS Level 3
209 */
210
211 /**
212 * Compute c <- a*b + beta * c (general matrix matrix
213 * multiplication)
214 */
215 public static DoubleMatrix gemm(double alpha, DoubleMatrix a,
216 DoubleMatrix b, double beta, DoubleMatrix c) {
217 Blas.dgemm('N', 'N', c.rows, c.columns, a.columns, alpha, a.data, 0,
218 a.rows, b.data, 0, b.rows, beta, c.data, 0, c.rows);
219 return c;
220 }
221
222 public static ComplexDoubleMatrix gemm(ComplexDouble alpha, ComplexDoubleMatrix a,
223 ComplexDoubleMatrix b, ComplexDouble beta, ComplexDoubleMatrix c) {
224 Blas.zgemm('N', 'N', c.rows, c.columns, a.columns, alpha, a.data, 0,
225 a.rows, b.data, 0, b.rows, beta, c.data, 0, c.rows);
226 return c;
227 }
228
229 /***************************************************************************
230 * LAPACK
231 */
232 public static DoubleMatrix gesv(DoubleMatrix a, int[] ipiv,
233 DoubleMatrix b) {
234 int info = Blas.dgesv(a.rows, b.columns, a.data, 0, a.rows, ipiv, 0,
235 b.data, 0, b.rows);
236 checkInfo("DGESV", info);
237
238 if (info > 0)
239 throw new LapackException("DGESV",
240 "Linear equation cannot be solved because the matrix was singular.");
241
242 return b;
243 }
244
245 //STOP
246 private static void checkInfo(String name, int info) {
247 if (info < -1)
248 throw new LapackArgumentException(name, info);
249 }
250 //START
251
252 public static DoubleMatrix sysv(char uplo, DoubleMatrix a, int[] ipiv,
253 DoubleMatrix b) {
254 int info = Blas.dsysv(uplo, a.rows, b.columns, a.data, 0, a.rows, ipiv, 0,
255 b.data, 0, b.rows);
256 checkInfo("DSYSV", info);
257
258 if (info > 0)
259 throw new IllegalArgumentException(
260 "Linear equation cannot be solved because the matrix was singular.");
261
262 return b;
263 }
264
265 public static int syev(char jobz, char uplo, DoubleMatrix a, DoubleMatrix w) {
266 double[] work = new double[1];
267 int info = Blas.dsyev(jobz, uplo, a.rows, a.data, 0, a.rows, w.data, 0,
268 work, 0, -1);
269 checkInfo("DSYEV", info);
270
271 int lwork = (int) work[0];
272 work = new double[lwork];
273
274 // System.out.println("Optimal LWORK = " + lwork);
275
276 info = Blas.dsyev(jobz, uplo, a.rows, a.data, 0, a.rows, w.data, 0,
277 work, 0, lwork);
278
279 if (info > 0)
280 throw new IllegalArgumentException(
281 "Eigenvalues could not be computed " + info
282 + " off-diagonal elements did not converge");
283
284 return info;
285 }
286
287 public static int syevx(char jobz, char range, char uplo, DoubleMatrix a,
288 double vl, double vu, int il, int iu, double abstol,
289 DoubleMatrix w, DoubleMatrix z) {
290 int n = a.rows;
291 int[] iwork = new int[5 * n];
292 int[] ifail = new int[n];
293 int[] m = new int[1];
294 int info;
295
296 info = Blas.dsyevx(jobz, range, uplo, n, a.data, 0, a.rows, vl, vu, il,
297 iu, abstol, m, 0, w.data, 0, z.data, 0, z.rows, iwork, 0, ifail, 0);
298
299 //System.out.printf("found eigenvalues = %d\n", m[0]);
300
301 if (info > 0) {
302 StringBuilder msg = new StringBuilder();
303 msg
304 .append("Not all eigenvalues converged. Non-converging eigenvalues were: ");
305 for (int i = 0; i < info; i++) {
306 if (i > 0)
307 msg.append(", ");
308 msg.append(ifail[i]);
309 }
310 msg.append(".");
311 throw new IllegalArgumentException(msg.toString());
312 }
313
314 return info;
315 }
316
317 public static int syevd(char jobz, char uplo, DoubleMatrix A,
318 DoubleMatrix w) {
319 int n = A.rows;
320
321 int info = Blas.dsyevd(jobz, uplo, n, A.data, 0, A.rows, w.data, 0);
322
323 if (info > 0)
324 throw new LapackConvergenceException("SYEVD", "Not all eigenvalues converged.");
325
326 return info;
327 }
328
329 public static int syevr(char jobz, char range, char uplo, DoubleMatrix a,
330 double vl, double vu, int il, int iu, double abstol,
331 DoubleMatrix w, DoubleMatrix z, int[] isuppz) {
332 int n = a.rows;
333 int[] m = new int[1];
334
335 int info = Blas.dsyevr(jobz, range, uplo, n, a.data, 0, a.rows, vl, vu,
336 il, iu, abstol, m, 0, w.data, 0, z.data, 0, z.rows, isuppz, 0);
337
338 //System.out.printf("found eigenvalues = %d\n", m[0]);
339
340 checkInfo("SYEVR", info);
341
342 return info;
343 }
344
345 public static void posv(char uplo, DoubleMatrix A, DoubleMatrix B) {
346 int n = A.rows;
347 int nrhs = B.columns;
348 int info = Blas.dposv(uplo, n, nrhs, A.data, 0, A.rows, B.data, 0,
349 B.rows);
350 checkInfo("DPOSV", info);
351 if (info > 0)
352 throw new LapackException("DPOSV",
353 "Leading minor of order i of A is not positive definite.");
354 }
355
356 public static int geev(char jobvl, char jobvr, DoubleMatrix A,
357 DoubleMatrix WR, DoubleMatrix WI, DoubleMatrix VL, DoubleMatrix VR) {
358 int info = Blas.dgeev(jobvl, jobvr, A.rows, A.data, 0, A.rows, WR.data, 0,
359 WI.data, 0, VL.data, 0, VL.rows, VR.data, 0, VR.rows);
360 if (info > 0)
361 throw new LapackException("DGEEV", "First " + info + " eigenvalues have not converged.");
362 return info;
363 }
364
365 //BEGIN
366 // The code below has been automatically generated.
367 // DO NOT EDIT!
368 /***************************************************************************
369 * BLAS Level 1
370 */
371
372 /** Compute x <-> y (swap two matrices) */
373 public static FloatMatrix swap(FloatMatrix x, FloatMatrix y) {
374 //Blas.sswap(x.length, x.data, 0, 1, y.data, 0, 1);
375 JavaBlas.rswap(x.length, x.data, 0, 1, y.data, 0, 1);
376 return y;
377 }
378
379 /** Compute x <- alpha * x (scale a matrix) */
380 public static FloatMatrix scal(float alpha, FloatMatrix x) {
381 Blas.sscal(x.length, alpha, x.data, 0, 1);
382 return x;
383 }
384
385 public static ComplexFloatMatrix scal(ComplexFloat alpha, ComplexFloatMatrix x) {
386 Blas.cscal(x.length, alpha, x.data, 0, 1);
387 return x;
388 }
389
390 /** Compute y <- x (copy a matrix) */
391 public static FloatMatrix copy(FloatMatrix x, FloatMatrix y) {
392 //Blas.scopy(x.length, x.data, 0, 1, y.data, 0, 1);
393 JavaBlas.rcopy(x.length, x.data, 0, 1, y.data, 0, 1);
394 return y;
395 }
396
397 public static ComplexFloatMatrix copy(ComplexFloatMatrix x, ComplexFloatMatrix y) {
398 Blas.ccopy(x.length, x.data, 0, 1, y.data, 0, 1);
399 return y;
400 }
401
402 /** Compute y <- alpha * x + y (elementwise addition) */
403 public static FloatMatrix axpy(float da, FloatMatrix dx, FloatMatrix dy) {
404 //Blas.saxpy(dx.length, da, dx.data, 0, 1, dy.data, 0, 1);
405 JavaBlas.raxpy(dx.length, da, dx.data, 0, 1, dy.data, 0, 1);
406
407 return dy;
408 }
409
410 public static ComplexFloatMatrix axpy(ComplexFloat da, ComplexFloatMatrix dx, ComplexFloatMatrix dy) {
411 Blas.caxpy(dx.length, da, dx.data, 0, 1, dy.data, 0, 1);
412 return dy;
413 }
414
415 /** Compute x^T * y (dot product) */
416 public static float dot(FloatMatrix x, FloatMatrix y) {
417 //return Blas.sdot(x.length, x.data, 0, 1, y.data, 0, 1);
418 return JavaBlas.rdot(x.length, x.data, 0, 1, y.data, 0, 1);
419 }
420
421 /** Compute x^T * y (dot product) */
422 public static ComplexFloat dotc(ComplexFloatMatrix x, ComplexFloatMatrix y) {
423 return Blas.cdotc(x.length, x.data, 0, 1, y.data, 0, 1);
424 }
425
426 /** Compute x^T * y (dot product) */
427 public static ComplexFloat dotu(ComplexFloatMatrix x, ComplexFloatMatrix y) {
428 return Blas.cdotu(x.length, x.data, 0, 1, y.data, 0, 1);
429 }
430
431 /** Compute || x ||_2 (2-norm) */
432 public static float nrm2(FloatMatrix x) {
433 return Blas.snrm2(x.length, x.data, 0, 1);
434 }
435
436 public static float nrm2(ComplexFloatMatrix x) {
437 return Blas.scnrm2(x.length, x.data, 0, 1);
438 }
439
440 /** Compute || x ||_1 (1-norm, sum of absolute values) */
441 public static float asum(FloatMatrix x) {
442 return Blas.sasum(x.length, x.data, 0, 1);
443 }
444
445 public static float asum(ComplexFloatMatrix x) {
446 return Blas.scasum(x.length, x.data, 0, 1);
447 }
448
449 /**
450 * Compute index of element with largest absolute value (index of absolute
451 * value maximum)
452 */
453 public static int iamax(FloatMatrix x) {
454 return Blas.isamax(x.length, x.data, 0, 1) - 1;
455 }
456
457 public static int iamax(ComplexFloatMatrix x) {
458 return Blas.icamax(x.length, x.data, 0, 1);
459 }
460
461 /***************************************************************************
462 * BLAS Level 2
463 */
464
465 /**
466 * Compute y <- alpha*op(a)*x + beta * y (general matrix vector
467 * multiplication)
468 */
469 public static FloatMatrix gemv(float alpha, FloatMatrix a,
470 FloatMatrix x, float beta, FloatMatrix y) {
471 if (false) {
472 Blas.sgemv('N', a.rows, a.columns, alpha, a.data, 0, a.rows, x.data, 0,
473 1, beta, y.data, 0, 1);
474 }
475 else {
476 if (beta == 0.0f) {
477 for (int i = 0; i < y.length; i++)
478 y.data[i] = 0.0f;
479
480 for (int j = 0; j < a.columns; j++)
481 for (int i = 0; i < a.rows; i++)
482 y.data[i] += a.get(i, j) * x.get(j);
483 }
484 else {
485 for (int j = 0; j < a.columns; j++)
486 for (int i = 0; i < a.rows; i++)
487 y.data[j] = a.get(i, j) * x.get(i) + y.data[j];
488 }
489 }
490 return y;
491 }
492
493 /** Compute A <- alpha * x * y^T + A (general rank-1 update) */
494 public static FloatMatrix ger(float alpha, FloatMatrix x,
495 FloatMatrix y, FloatMatrix a) {
496 Blas.sger(a.rows, a.columns, alpha, x.data, 0, 1, y.data, 0, 1, a.data,
497 0, a.rows);
498 return a;
499 }
500
501 /** Compute A <- alpha * x * y^T + A (general rank-1 update) */
502 public static ComplexFloatMatrix geru(ComplexFloat alpha, ComplexFloatMatrix x,
503 ComplexFloatMatrix y, ComplexFloatMatrix a) {
504 Blas.cgeru(a.rows, a.columns, alpha, x.data, 0, 1, y.data, 0, 1, a.data,
505 0, a.rows);
506 return a;
507 }
508
509 /** Compute A <- alpha * x * y^H + A (general rank-1 update) */
510 public static ComplexFloatMatrix gerc(ComplexFloat alpha, ComplexFloatMatrix x,
511 ComplexFloatMatrix y, ComplexFloatMatrix a) {
512 Blas.cgerc(a.rows, a.columns, alpha, x.data, 0, 1, y.data, 0, 1, a.data,
513 0, a.rows);
514 return a;
515 }
516
517 /***************************************************************************
518 * BLAS Level 3
519 */
520
521 /**
522 * Compute c <- a*b + beta * c (general matrix matrix
523 * multiplication)
524 */
525 public static FloatMatrix gemm(float alpha, FloatMatrix a,
526 FloatMatrix b, float beta, FloatMatrix c) {
527 Blas.sgemm('N', 'N', c.rows, c.columns, a.columns, alpha, a.data, 0,
528 a.rows, b.data, 0, b.rows, beta, c.data, 0, c.rows);
529 return c;
530 }
531
532 public static ComplexFloatMatrix gemm(ComplexFloat alpha, ComplexFloatMatrix a,
533 ComplexFloatMatrix b, ComplexFloat beta, ComplexFloatMatrix c) {
534 Blas.cgemm('N', 'N', c.rows, c.columns, a.columns, alpha, a.data, 0,
535 a.rows, b.data, 0, b.rows, beta, c.data, 0, c.rows);
536 return c;
537 }
538
539 /***************************************************************************
540 * LAPACK
541 */
542 public static FloatMatrix gesv(FloatMatrix a, int[] ipiv,
543 FloatMatrix b) {
544 int info = Blas.sgesv(a.rows, b.columns, a.data, 0, a.rows, ipiv, 0,
545 b.data, 0, b.rows);
546 checkInfo("DGESV", info);
547
548 if (info > 0)
549 throw new LapackException("DGESV",
550 "Linear equation cannot be solved because the matrix was singular.");
551
552 return b;
553 }
554
555
556 public static FloatMatrix sysv(char uplo, FloatMatrix a, int[] ipiv,
557 FloatMatrix b) {
558 int info = Blas.ssysv(uplo, a.rows, b.columns, a.data, 0, a.rows, ipiv, 0,
559 b.data, 0, b.rows);
560 checkInfo("DSYSV", info);
561
562 if (info > 0)
563 throw new IllegalArgumentException(
564 "Linear equation cannot be solved because the matrix was singular.");
565
566 return b;
567 }
568
569 public static int syev(char jobz, char uplo, FloatMatrix a, FloatMatrix w) {
570 float[] work = new float[1];
571 int info = Blas.ssyev(jobz, uplo, a.rows, a.data, 0, a.rows, w.data, 0,
572 work, 0, -1);
573 checkInfo("DSYEV", info);
574
575 int lwork = (int) work[0];
576 work = new float[lwork];
577
578 // System.out.println("Optimal LWORK = " + lwork);
579
580 info = Blas.ssyev(jobz, uplo, a.rows, a.data, 0, a.rows, w.data, 0,
581 work, 0, lwork);
582
583 if (info > 0)
584 throw new IllegalArgumentException(
585 "Eigenvalues could not be computed " + info
586 + " off-diagonal elements did not converge");
587
588 return info;
589 }
590
591 public static int syevx(char jobz, char range, char uplo, FloatMatrix a,
592 float vl, float vu, int il, int iu, float abstol,
593 FloatMatrix w, FloatMatrix z) {
594 int n = a.rows;
595 int[] iwork = new int[5 * n];
596 int[] ifail = new int[n];
597 int[] m = new int[1];
598 int info;
599
600 info = Blas.ssyevx(jobz, range, uplo, n, a.data, 0, a.rows, vl, vu, il,
601 iu, abstol, m, 0, w.data, 0, z.data, 0, z.rows, iwork, 0, ifail, 0);
602
603 //System.out.printf("found eigenvalues = %d\n", m[0]);
604
605 if (info > 0) {
606 StringBuilder msg = new StringBuilder();
607 msg
608 .append("Not all eigenvalues converged. Non-converging eigenvalues were: ");
609 for (int i = 0; i < info; i++) {
610 if (i > 0)
611 msg.append(", ");
612 msg.append(ifail[i]);
613 }
614 msg.append(".");
615 throw new IllegalArgumentException(msg.toString());
616 }
617
618 return info;
619 }
620
621 public static int syevd(char jobz, char uplo, FloatMatrix A,
622 FloatMatrix w) {
623 int n = A.rows;
624
625 int info = Blas.ssyevd(jobz, uplo, n, A.data, 0, A.rows, w.data, 0);
626
627 if (info > 0)
628 throw new LapackConvergenceException("SYEVD", "Not all eigenvalues converged.");
629
630 return info;
631 }
632
633 public static int syevr(char jobz, char range, char uplo, FloatMatrix a,
634 float vl, float vu, int il, int iu, float abstol,
635 FloatMatrix w, FloatMatrix z, int[] isuppz) {
636 int n = a.rows;
637 int[] m = new int[1];
638
639 int info = Blas.ssyevr(jobz, range, uplo, n, a.data, 0, a.rows, vl, vu,
640 il, iu, abstol, m, 0, w.data, 0, z.data, 0, z.rows, isuppz, 0);
641
642 //System.out.printf("found eigenvalues = %d\n", m[0]);
643
644 checkInfo("SYEVR", info);
645
646 return info;
647 }
648
649 public static void posv(char uplo, FloatMatrix A, FloatMatrix B) {
650 int n = A.rows;
651 int nrhs = B.columns;
652 int info = Blas.sposv(uplo, n, nrhs, A.data, 0, A.rows, B.data, 0,
653 B.rows);
654 checkInfo("DPOSV", info);
655 if (info > 0)
656 throw new LapackException("DPOSV",
657 "Leading minor of order i of A is not positive definite.");
658 }
659
660 public static int geev(char jobvl, char jobvr, FloatMatrix A,
661 FloatMatrix WR, FloatMatrix WI, FloatMatrix VL, FloatMatrix VR) {
662 int info = Blas.sgeev(jobvl, jobvr, A.rows, A.data, 0, A.rows, WR.data, 0,
663 WI.data, 0, VL.data, 0, VL.rows, VR.data, 0, VR.rows);
664 if (info > 0)
665 throw new LapackException("DGEEV", "First " + info + " eigenvalues have not converged.");
666 return info;
667 }
668
669 //END
670 }