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 &lt;-&gt; 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 &lt;-&gt; 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}