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