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 }