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 }