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
038package org.jblas;
039
040import org.jblas.exceptions.LapackException;
041import org.jblas.exceptions.LapackArgumentException;
042import org.jblas.exceptions.LapackConvergenceException;
043import 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 */
056public 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}