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