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 }