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