001    // --- BEGIN LICENSE BLOCK ---
002    /* 
003     * Copyright (c) 2009-2011, Mikio L. Braun
004     *               2011, Nicolas Oury
005     * All rights reserved.
006     * 
007     * Redistribution and use in source and binary forms, with or without
008     * modification, are permitted provided that the following conditions are
009     * met:
010     * 
011     *     * Redistributions of source code must retain the above copyright
012     *       notice, this list of conditions and the following disclaimer.
013     * 
014     *     * Redistributions in binary form must reproduce the above
015     *       copyright notice, this list of conditions and the following
016     *       disclaimer in the documentation and/or other materials provided
017     *       with the distribution.
018     * 
019     *     * Neither the name of the Technische Universit?t Berlin nor the
020     *       names of its contributors may be used to endorse or promote
021     *       products derived from this software without specific prior
022     *       written permission.
023     * 
024     * THIS SOFTWARE IS PROVIDED BY THE COPYRIGHT HOLDERS AND CONTRIBUTORS
025     * "AS IS" AND ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT
026     * LIMITED TO, THE IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR
027     * A PARTICULAR PURPOSE ARE DISCLAIMED. IN NO EVENT SHALL THE COPYRIGHT
028     * HOLDER OR CONTRIBUTORS BE LIABLE FOR ANY DIRECT, INDIRECT, INCIDENTAL,
029     * SPECIAL, EXEMPLARY, OR CONSEQUENTIAL DAMAGES (INCLUDING, BUT NOT
030     * LIMITED TO, PROCUREMENT OF SUBSTITUTE GOODS OR SERVICES; LOSS OF USE,
031     * DATA, OR PROFITS; OR BUSINESS INTERRUPTION) HOWEVER CAUSED AND ON ANY
032     * THEORY OF LIABILITY, WHETHER IN CONTRACT, STRICT LIABILITY, OR TORT
033     * (INCLUDING NEGLIGENCE OR OTHERWISE) ARISING IN ANY WAY OUT OF THE USE
034     * OF THIS SOFTWARE, EVEN IF ADVISED OF THE POSSIBILITY OF SUCH DAMAGE.
035     */
036    // --- END LICENSE BLOCK ---
037    
038    package org.jblas;
039    
040    import org.jblas.exceptions.LapackException;
041    import org.jblas.exceptions.LapackArgumentException;
042    import org.jblas.exceptions.LapackConvergenceException;
043    import org.jblas.exceptions.LapackSingularityException;
044    
045    //import edu.ida.core.OutputValue;
046    
047    /**
048     * This class provides a cleaner direct interface to the BLAS routines by
049     * extracting the parameters of the matrices from the matrices itself.
050     * <p/>
051     * For example, you can just pass the vector and do not have to pass the length,
052     * corresponding DoubleBuffer, offset and step size explicitly.
053     * <p/>
054     * Currently, all the general matrix routines are implemented.
055     */
056    public class SimpleBlas {
057        /***************************************************************************
058         * BLAS Level 1
059         */
060    
061        /**
062         * Compute x <-> y (swap two matrices)
063         */
064        public static DoubleMatrix swap(DoubleMatrix x, DoubleMatrix y) {
065            //NativeBlas.dswap(x.length, x.data, 0, 1, y.data, 0, 1);
066            JavaBlas.rswap(x.length, x.data, 0, 1, y.data, 0, 1);
067            return y;
068        }
069    
070        /**
071         * Compute x <- alpha * x (scale a matrix)
072         */
073        public static DoubleMatrix scal(double alpha, DoubleMatrix x) {
074            NativeBlas.dscal(x.length, alpha, x.data, 0, 1);
075            return x;
076        }
077    
078        public static ComplexDoubleMatrix scal(ComplexDouble alpha, ComplexDoubleMatrix x) {
079            NativeBlas.zscal(x.length, alpha, x.data, 0, 1);
080            return x;
081        }
082    
083        /**
084         * Compute y <- x (copy a matrix)
085         */
086        public static DoubleMatrix copy(DoubleMatrix x, DoubleMatrix y) {
087            //NativeBlas.dcopy(x.length, x.data, 0, 1, y.data, 0, 1);
088            JavaBlas.rcopy(x.length, x.data, 0, 1, y.data, 0, 1);
089            return y;
090        }
091    
092        public static ComplexDoubleMatrix copy(ComplexDoubleMatrix x, ComplexDoubleMatrix y) {
093            NativeBlas.zcopy(x.length, x.data, 0, 1, y.data, 0, 1);
094            return y;
095        }
096    
097        /**
098         * Compute y <- alpha * x + y (elementwise addition)
099         */
100        public static DoubleMatrix axpy(double da, DoubleMatrix dx, DoubleMatrix dy) {
101            //NativeBlas.daxpy(dx.length, da, dx.data, 0, 1, dy.data, 0, 1);
102            JavaBlas.raxpy(dx.length, da, dx.data, 0, 1, dy.data, 0, 1);
103    
104            return dy;
105        }
106    
107        public static ComplexDoubleMatrix axpy(ComplexDouble da, ComplexDoubleMatrix dx, ComplexDoubleMatrix dy) {
108            NativeBlas.zaxpy(dx.length, da, dx.data, 0, 1, dy.data, 0, 1);
109            return dy;
110        }
111    
112        /**
113         * Compute x^T * y (dot product)
114         */
115        public static double dot(DoubleMatrix x, DoubleMatrix y) {
116            //return NativeBlas.ddot(x.length, x.data, 0, 1, y.data, 0, 1);
117            return JavaBlas.rdot(x.length, x.data, 0, 1, y.data, 0, 1);
118        }
119    
120        /**
121         * Compute x^T * y (dot product)
122         */
123        public static ComplexDouble dotc(ComplexDoubleMatrix x, ComplexDoubleMatrix y) {
124            return NativeBlas.zdotc(x.length, x.data, 0, 1, y.data, 0, 1);
125        }
126    
127        /**
128         * Compute x^T * y (dot product)
129         */
130        public static ComplexDouble dotu(ComplexDoubleMatrix x, ComplexDoubleMatrix y) {
131            return NativeBlas.zdotu(x.length, x.data, 0, 1, y.data, 0, 1);
132        }
133    
134        /**
135         * Compute || x ||_2 (2-norm)
136         */
137        public static double nrm2(DoubleMatrix x) {
138            return NativeBlas.dnrm2(x.length, x.data, 0, 1);
139        }
140    
141        public static double nrm2(ComplexDoubleMatrix x) {
142            return NativeBlas.dznrm2(x.length, x.data, 0, 1);
143        }
144    
145        /**
146         * Compute || x ||_1 (1-norm, sum of absolute values)
147         */
148        public static double asum(DoubleMatrix x) {
149            return NativeBlas.dasum(x.length, x.data, 0, 1);
150        }
151    
152        public static double asum(ComplexDoubleMatrix x) {
153            return NativeBlas.dzasum(x.length, x.data, 0, 1);
154        }
155    
156        /**
157         * Compute index of element with largest absolute value (index of absolute
158         * value maximum)
159         */
160        public static int iamax(DoubleMatrix x) {
161            return NativeBlas.idamax(x.length, x.data, 0, 1) - 1;
162        }
163    
164        public static int iamax(ComplexDoubleMatrix x) {
165            return NativeBlas.izamax(x.length, x.data, 0, 1);
166        }
167    
168        /***************************************************************************
169         * BLAS Level 2
170         */
171    
172        /**
173         * Compute y <- alpha*op(a)*x + beta * y (general matrix vector
174         * multiplication)
175         */
176        public static DoubleMatrix gemv(double alpha, DoubleMatrix a,
177                                        DoubleMatrix x, double beta, DoubleMatrix y) {
178            if (false) {
179                NativeBlas.dgemv('N', a.rows, a.columns, alpha, a.data, 0, a.rows, x.data, 0,
180                        1, beta, y.data, 0, 1);
181            } else {
182                if (beta == 0.0) {
183                    for (int i = 0; i < y.length; i++)
184                        y.data[i] = 0.0;
185    
186                    for (int j = 0; j < a.columns; j++) {
187                        double xj = x.get(j);
188                        if (xj != 0.0) {
189                            for (int i = 0; i < a.rows; i++)
190                                y.data[i] += a.get(i, j) * xj;
191                        }
192                    }
193                } else {
194                    for (int j = 0; j < a.columns; j++) {
195                        double byj = beta * y.data[j];
196                        double xj = x.get(j);
197                        for (int i = 0; i < a.rows; i++)
198                            y.data[j] = a.get(i, j) * xj + byj;
199                    }
200                }
201            }
202            return y;
203        }
204    
205        /**
206         * Compute A <- alpha * x * y^T + A (general rank-1 update)
207         */
208        public static DoubleMatrix ger(double alpha, DoubleMatrix x,
209                                       DoubleMatrix y, DoubleMatrix a) {
210            NativeBlas.dger(a.rows, a.columns, alpha, x.data, 0, 1, y.data, 0, 1, a.data,
211                    0, a.rows);
212            return a;
213        }
214    
215        /**
216         * Compute A <- alpha * x * y^T + A (general rank-1 update)
217         */
218        public static ComplexDoubleMatrix geru(ComplexDouble alpha, ComplexDoubleMatrix x,
219                                               ComplexDoubleMatrix y, ComplexDoubleMatrix a) {
220            NativeBlas.zgeru(a.rows, a.columns, alpha, x.data, 0, 1, y.data, 0, 1, a.data,
221                    0, a.rows);
222            return a;
223        }
224    
225        /**
226         * Compute A <- alpha * x * y^H + A (general rank-1 update)
227         */
228        public static ComplexDoubleMatrix gerc(ComplexDouble alpha, ComplexDoubleMatrix x,
229                                               ComplexDoubleMatrix y, ComplexDoubleMatrix a) {
230            NativeBlas.zgerc(a.rows, a.columns, alpha, x.data, 0, 1, y.data, 0, 1, a.data,
231                    0, a.rows);
232            return a;
233        }
234    
235        /***************************************************************************
236         * BLAS Level 3
237         */
238    
239        /**
240         * Compute c <- a*b + beta * c (general matrix matrix
241         * multiplication)
242         */
243        public static DoubleMatrix gemm(double alpha, DoubleMatrix a,
244                                        DoubleMatrix b, double beta, DoubleMatrix c) {
245            NativeBlas.dgemm('N', 'N', c.rows, c.columns, a.columns, alpha, a.data, 0,
246                    a.rows, b.data, 0, b.rows, beta, c.data, 0, c.rows);
247            return c;
248        }
249    
250        public static ComplexDoubleMatrix gemm(ComplexDouble alpha, ComplexDoubleMatrix a,
251                                               ComplexDoubleMatrix b, ComplexDouble beta, ComplexDoubleMatrix c) {
252            NativeBlas.zgemm('N', 'N', c.rows, c.columns, a.columns, alpha, a.data, 0,
253                    a.rows, b.data, 0, b.rows, beta, c.data, 0, c.rows);
254            return c;
255        }
256    
257            /***************************************************************************
258         * LAPACK
259         */
260    
261        public static DoubleMatrix gesv(DoubleMatrix a, int[] ipiv,
262                                        DoubleMatrix b) {
263            int info = NativeBlas.dgesv(a.rows, b.columns, a.data, 0, a.rows, ipiv, 0,
264                    b.data, 0, b.rows);
265            checkInfo("DGESV", info);
266    
267            if (info > 0)
268                throw new LapackException("DGESV",
269                        "Linear equation cannot be solved because the matrix was singular.");
270    
271            return b;
272        }
273    
274    //STOP
275    
276        private static void checkInfo(String name, int info) {
277            if (info < -1)
278                throw new LapackArgumentException(name, info);
279        }
280    //START
281    
282        public static DoubleMatrix sysv(char uplo, DoubleMatrix a, int[] ipiv,
283                                        DoubleMatrix b) {
284            int info = NativeBlas.dsysv(uplo, a.rows, b.columns, a.data, 0, a.rows, ipiv, 0,
285                    b.data, 0, b.rows);
286            checkInfo("SYSV", info);
287    
288            if (info > 0)
289                throw new LapackSingularityException("SYV",
290                        "Linear equation cannot be solved because the matrix was singular.");
291    
292            return b;
293        }
294    
295        public static int syev(char jobz, char uplo, DoubleMatrix a, DoubleMatrix w) {
296            int info = NativeBlas.dsyev(jobz, uplo, a.rows, a.data, 0, a.rows, w.data, 0);
297    
298            if (info > 0)
299                throw new LapackConvergenceException("SYEV",
300                        "Eigenvalues could not be computed " + info
301                                + " off-diagonal elements did not converge");
302    
303            return info;
304        }
305    
306        public static int syevx(char jobz, char range, char uplo, DoubleMatrix a,
307                                double vl, double vu, int il, int iu, double abstol,
308                                DoubleMatrix w, DoubleMatrix z) {
309            int n = a.rows;
310            int[] iwork = new int[5 * n];
311            int[] ifail = new int[n];
312            int[] m = new int[1];
313            int info;
314    
315            info = NativeBlas.dsyevx(jobz, range, uplo, n, a.data, 0, a.rows, vl, vu, il,
316                    iu, abstol, m, 0, w.data, 0, z.data, 0, z.rows, iwork, 0, ifail, 0);
317    
318            if (info > 0) {
319                StringBuilder msg = new StringBuilder();
320                msg
321                        .append("Not all eigenvalues converged. Non-converging eigenvalues were: ");
322                for (int i = 0; i < info; i++) {
323                    if (i > 0)
324                        msg.append(", ");
325                    msg.append(ifail[i]);
326                }
327                msg.append(".");
328                throw new LapackConvergenceException("SYEVX", msg.toString());
329            }
330    
331            return info;
332        }
333    
334        public static int syevd(char jobz, char uplo, DoubleMatrix A,
335                                DoubleMatrix w) {
336            int n = A.rows;
337    
338            int info = NativeBlas.dsyevd(jobz, uplo, n, A.data, 0, A.rows, w.data, 0);
339    
340            if (info > 0)
341                throw new LapackConvergenceException("SYEVD", "Not all eigenvalues converged.");
342    
343            return info;
344        }
345    
346        public static int syevr(char jobz, char range, char uplo, DoubleMatrix a,
347                                double vl, double vu, int il, int iu, double abstol,
348                                DoubleMatrix w, DoubleMatrix z, int[] isuppz) {
349            int n = a.rows;
350            int[] m = new int[1];
351    
352            int info = NativeBlas.dsyevr(jobz, range, uplo, n, a.data, 0, a.rows, vl, vu,
353                    il, iu, abstol, m, 0, w.data, 0, z.data, 0, z.rows, isuppz, 0);
354    
355            checkInfo("SYEVR", info);
356    
357            return info;
358        }
359    
360        public static void posv(char uplo, DoubleMatrix A, DoubleMatrix B) {
361            int n = A.rows;
362            int nrhs = B.columns;
363            int info = NativeBlas.dposv(uplo, n, nrhs, A.data, 0, A.rows, B.data, 0,
364                    B.rows);
365            checkInfo("DPOSV", info);
366            if (info > 0)
367                throw new LapackArgumentException("DPOSV",
368                        "Leading minor of order i of A is not positive definite.");
369        }
370    
371        public static int geev(char jobvl, char jobvr, DoubleMatrix A,
372                               DoubleMatrix WR, DoubleMatrix WI, DoubleMatrix VL, DoubleMatrix VR) {
373            int info = NativeBlas.dgeev(jobvl, jobvr, A.rows, A.data, 0, A.rows, WR.data, 0,
374                    WI.data, 0, VL.data, 0, VL.rows, VR.data, 0, VR.rows);
375            if (info > 0)
376                throw new LapackConvergenceException("DGEEV", "First " + info + " eigenvalues have not converged.");
377            return info;
378        }
379    
380        public static int sygvd(int itype, char jobz, char uplo, DoubleMatrix A, DoubleMatrix B, DoubleMatrix W) {
381            int info = NativeBlas.dsygvd(itype, jobz, uplo, A.rows, A.data, 0, A.rows, B.data, 0, B.rows, W.data, 0);
382            if (info == 0)
383                return 0;
384            else {
385                if (info < 0)
386                    throw new LapackArgumentException("DSYGVD", -info);
387                if (info <= A.rows && jobz == 'N')
388                    throw new LapackConvergenceException("DSYGVD", info + " off-diagonal elements did not converge to 0.");
389                if (info <= A.rows && jobz == 'V')
390                    throw new LapackException("DSYGVD", "Failed to compute an eigenvalue while working on a sub-matrix  " + info + ".");
391                else
392                    throw new LapackException("DSYGVD", "The leading minor of order " + (info - A.rows) + " of B is not positive definite.");
393            }
394        }
395    
396    //BEGIN
397      // The code below has been automatically generated.
398      // DO NOT EDIT!
399        /***************************************************************************
400         * BLAS Level 1
401         */
402    
403        /**
404         * Compute x <-> y (swap two matrices)
405         */
406        public static FloatMatrix swap(FloatMatrix x, FloatMatrix y) {
407            //NativeBlas.sswap(x.length, x.data, 0, 1, y.data, 0, 1);
408            JavaBlas.rswap(x.length, x.data, 0, 1, y.data, 0, 1);
409            return y;
410        }
411    
412        /**
413         * Compute x <- alpha * x (scale a matrix)
414         */
415        public static FloatMatrix scal(float alpha, FloatMatrix x) {
416            NativeBlas.sscal(x.length, alpha, x.data, 0, 1);
417            return x;
418        }
419    
420        public static ComplexFloatMatrix scal(ComplexFloat alpha, ComplexFloatMatrix x) {
421            NativeBlas.cscal(x.length, alpha, x.data, 0, 1);
422            return x;
423        }
424    
425        /**
426         * Compute y <- x (copy a matrix)
427         */
428        public static FloatMatrix copy(FloatMatrix x, FloatMatrix y) {
429            //NativeBlas.scopy(x.length, x.data, 0, 1, y.data, 0, 1);
430            JavaBlas.rcopy(x.length, x.data, 0, 1, y.data, 0, 1);
431            return y;
432        }
433    
434        public static ComplexFloatMatrix copy(ComplexFloatMatrix x, ComplexFloatMatrix y) {
435            NativeBlas.ccopy(x.length, x.data, 0, 1, y.data, 0, 1);
436            return y;
437        }
438    
439        /**
440         * Compute y <- alpha * x + y (elementwise addition)
441         */
442        public static FloatMatrix axpy(float da, FloatMatrix dx, FloatMatrix dy) {
443            //NativeBlas.saxpy(dx.length, da, dx.data, 0, 1, dy.data, 0, 1);
444            JavaBlas.raxpy(dx.length, da, dx.data, 0, 1, dy.data, 0, 1);
445    
446            return dy;
447        }
448    
449        public static ComplexFloatMatrix axpy(ComplexFloat da, ComplexFloatMatrix dx, ComplexFloatMatrix dy) {
450            NativeBlas.caxpy(dx.length, da, dx.data, 0, 1, dy.data, 0, 1);
451            return dy;
452        }
453    
454        /**
455         * Compute x^T * y (dot product)
456         */
457        public static float dot(FloatMatrix x, FloatMatrix y) {
458            //return NativeBlas.sdot(x.length, x.data, 0, 1, y.data, 0, 1);
459            return JavaBlas.rdot(x.length, x.data, 0, 1, y.data, 0, 1);
460        }
461    
462        /**
463         * Compute x^T * y (dot product)
464         */
465        public static ComplexFloat dotc(ComplexFloatMatrix x, ComplexFloatMatrix y) {
466            return NativeBlas.cdotc(x.length, x.data, 0, 1, y.data, 0, 1);
467        }
468    
469        /**
470         * Compute x^T * y (dot product)
471         */
472        public static ComplexFloat dotu(ComplexFloatMatrix x, ComplexFloatMatrix y) {
473            return NativeBlas.cdotu(x.length, x.data, 0, 1, y.data, 0, 1);
474        }
475    
476        /**
477         * Compute || x ||_2 (2-norm)
478         */
479        public static float nrm2(FloatMatrix x) {
480            return NativeBlas.snrm2(x.length, x.data, 0, 1);
481        }
482    
483        public static float nrm2(ComplexFloatMatrix x) {
484            return NativeBlas.scnrm2(x.length, x.data, 0, 1);
485        }
486    
487        /**
488         * Compute || x ||_1 (1-norm, sum of absolute values)
489         */
490        public static float asum(FloatMatrix x) {
491            return NativeBlas.sasum(x.length, x.data, 0, 1);
492        }
493    
494        public static float asum(ComplexFloatMatrix x) {
495            return NativeBlas.scasum(x.length, x.data, 0, 1);
496        }
497    
498        /**
499         * Compute index of element with largest absolute value (index of absolute
500         * value maximum)
501         */
502        public static int iamax(FloatMatrix x) {
503            return NativeBlas.isamax(x.length, x.data, 0, 1) - 1;
504        }
505    
506        public static int iamax(ComplexFloatMatrix x) {
507            return NativeBlas.icamax(x.length, x.data, 0, 1);
508        }
509    
510        /***************************************************************************
511         * BLAS Level 2
512         */
513    
514        /**
515         * Compute y <- alpha*op(a)*x + beta * y (general matrix vector
516         * multiplication)
517         */
518        public static FloatMatrix gemv(float alpha, FloatMatrix a,
519                                        FloatMatrix x, float beta, FloatMatrix y) {
520            if (false) {
521                NativeBlas.sgemv('N', a.rows, a.columns, alpha, a.data, 0, a.rows, x.data, 0,
522                        1, beta, y.data, 0, 1);
523            } else {
524                if (beta == 0.0f) {
525                    for (int i = 0; i < y.length; i++)
526                        y.data[i] = 0.0f;
527    
528                    for (int j = 0; j < a.columns; j++) {
529                        float xj = x.get(j);
530                        if (xj != 0.0f) {
531                            for (int i = 0; i < a.rows; i++)
532                                y.data[i] += a.get(i, j) * xj;
533                        }
534                    }
535                } else {
536                    for (int j = 0; j < a.columns; j++) {
537                        float byj = beta * y.data[j];
538                        float xj = x.get(j);
539                        for (int i = 0; i < a.rows; i++)
540                            y.data[j] = a.get(i, j) * xj + byj;
541                    }
542                }
543            }
544            return y;
545        }
546    
547        /**
548         * Compute A <- alpha * x * y^T + A (general rank-1 update)
549         */
550        public static FloatMatrix ger(float alpha, FloatMatrix x,
551                                       FloatMatrix y, FloatMatrix a) {
552            NativeBlas.sger(a.rows, a.columns, alpha, x.data, 0, 1, y.data, 0, 1, a.data,
553                    0, a.rows);
554            return a;
555        }
556    
557        /**
558         * Compute A <- alpha * x * y^T + A (general rank-1 update)
559         */
560        public static ComplexFloatMatrix geru(ComplexFloat alpha, ComplexFloatMatrix x,
561                                               ComplexFloatMatrix y, ComplexFloatMatrix a) {
562            NativeBlas.cgeru(a.rows, a.columns, alpha, x.data, 0, 1, y.data, 0, 1, a.data,
563                    0, a.rows);
564            return a;
565        }
566    
567        /**
568         * Compute A <- alpha * x * y^H + A (general rank-1 update)
569         */
570        public static ComplexFloatMatrix gerc(ComplexFloat alpha, ComplexFloatMatrix x,
571                                               ComplexFloatMatrix y, ComplexFloatMatrix a) {
572            NativeBlas.cgerc(a.rows, a.columns, alpha, x.data, 0, 1, y.data, 0, 1, a.data,
573                    0, a.rows);
574            return a;
575        }
576    
577        /***************************************************************************
578         * BLAS Level 3
579         */
580    
581        /**
582         * Compute c <- a*b + beta * c (general matrix matrix
583         * multiplication)
584         */
585        public static FloatMatrix gemm(float alpha, FloatMatrix a,
586                                        FloatMatrix b, float beta, FloatMatrix c) {
587            NativeBlas.sgemm('N', 'N', c.rows, c.columns, a.columns, alpha, a.data, 0,
588                    a.rows, b.data, 0, b.rows, beta, c.data, 0, c.rows);
589            return c;
590        }
591    
592        public static ComplexFloatMatrix gemm(ComplexFloat alpha, ComplexFloatMatrix a,
593                                               ComplexFloatMatrix b, ComplexFloat beta, ComplexFloatMatrix c) {
594            NativeBlas.cgemm('N', 'N', c.rows, c.columns, a.columns, alpha, a.data, 0,
595                    a.rows, b.data, 0, b.rows, beta, c.data, 0, c.rows);
596            return c;
597        }
598    
599            /***************************************************************************
600         * LAPACK
601         */
602    
603        public static FloatMatrix gesv(FloatMatrix a, int[] ipiv,
604                                        FloatMatrix b) {
605            int info = NativeBlas.sgesv(a.rows, b.columns, a.data, 0, a.rows, ipiv, 0,
606                    b.data, 0, b.rows);
607            checkInfo("DGESV", info);
608    
609            if (info > 0)
610                throw new LapackException("DGESV",
611                        "Linear equation cannot be solved because the matrix was singular.");
612    
613            return b;
614        }
615    
616    
617        public static FloatMatrix sysv(char uplo, FloatMatrix a, int[] ipiv,
618                                        FloatMatrix b) {
619            int info = NativeBlas.ssysv(uplo, a.rows, b.columns, a.data, 0, a.rows, ipiv, 0,
620                    b.data, 0, b.rows);
621            checkInfo("SYSV", info);
622    
623            if (info > 0)
624                throw new LapackSingularityException("SYV",
625                        "Linear equation cannot be solved because the matrix was singular.");
626    
627            return b;
628        }
629    
630        public static int syev(char jobz, char uplo, FloatMatrix a, FloatMatrix w) {
631            int info = NativeBlas.ssyev(jobz, uplo, a.rows, a.data, 0, a.rows, w.data, 0);
632    
633            if (info > 0)
634                throw new LapackConvergenceException("SYEV",
635                        "Eigenvalues could not be computed " + info
636                                + " off-diagonal elements did not converge");
637    
638            return info;
639        }
640    
641        public static int syevx(char jobz, char range, char uplo, FloatMatrix a,
642                                float vl, float vu, int il, int iu, float abstol,
643                                FloatMatrix w, FloatMatrix z) {
644            int n = a.rows;
645            int[] iwork = new int[5 * n];
646            int[] ifail = new int[n];
647            int[] m = new int[1];
648            int info;
649    
650            info = NativeBlas.ssyevx(jobz, range, uplo, n, a.data, 0, a.rows, vl, vu, il,
651                    iu, abstol, m, 0, w.data, 0, z.data, 0, z.rows, iwork, 0, ifail, 0);
652    
653            if (info > 0) {
654                StringBuilder msg = new StringBuilder();
655                msg
656                        .append("Not all eigenvalues converged. Non-converging eigenvalues were: ");
657                for (int i = 0; i < info; i++) {
658                    if (i > 0)
659                        msg.append(", ");
660                    msg.append(ifail[i]);
661                }
662                msg.append(".");
663                throw new LapackConvergenceException("SYEVX", msg.toString());
664            }
665    
666            return info;
667        }
668    
669        public static int syevd(char jobz, char uplo, FloatMatrix A,
670                                FloatMatrix w) {
671            int n = A.rows;
672    
673            int info = NativeBlas.ssyevd(jobz, uplo, n, A.data, 0, A.rows, w.data, 0);
674    
675            if (info > 0)
676                throw new LapackConvergenceException("SYEVD", "Not all eigenvalues converged.");
677    
678            return info;
679        }
680    
681        public static int syevr(char jobz, char range, char uplo, FloatMatrix a,
682                                float vl, float vu, int il, int iu, float abstol,
683                                FloatMatrix w, FloatMatrix z, int[] isuppz) {
684            int n = a.rows;
685            int[] m = new int[1];
686    
687            int info = NativeBlas.ssyevr(jobz, range, uplo, n, a.data, 0, a.rows, vl, vu,
688                    il, iu, abstol, m, 0, w.data, 0, z.data, 0, z.rows, isuppz, 0);
689    
690            checkInfo("SYEVR", info);
691    
692            return info;
693        }
694    
695        public static void posv(char uplo, FloatMatrix A, FloatMatrix B) {
696            int n = A.rows;
697            int nrhs = B.columns;
698            int info = NativeBlas.sposv(uplo, n, nrhs, A.data, 0, A.rows, B.data, 0,
699                    B.rows);
700            checkInfo("DPOSV", info);
701            if (info > 0)
702                throw new LapackArgumentException("DPOSV",
703                        "Leading minor of order i of A is not positive definite.");
704        }
705    
706        public static int geev(char jobvl, char jobvr, FloatMatrix A,
707                               FloatMatrix WR, FloatMatrix WI, FloatMatrix VL, FloatMatrix VR) {
708            int info = NativeBlas.sgeev(jobvl, jobvr, A.rows, A.data, 0, A.rows, WR.data, 0,
709                    WI.data, 0, VL.data, 0, VL.rows, VR.data, 0, VR.rows);
710            if (info > 0)
711                throw new LapackConvergenceException("DGEEV", "First " + info + " eigenvalues have not converged.");
712            return info;
713        }
714    
715        public static int sygvd(int itype, char jobz, char uplo, FloatMatrix A, FloatMatrix B, FloatMatrix W) {
716            int info = NativeBlas.ssygvd(itype, jobz, uplo, A.rows, A.data, 0, A.rows, B.data, 0, B.rows, W.data, 0);
717            if (info == 0)
718                return 0;
719            else {
720                if (info < 0)
721                    throw new LapackArgumentException("DSYGVD", -info);
722                if (info <= A.rows && jobz == 'N')
723                    throw new LapackConvergenceException("DSYGVD", info + " off-diagonal elements did not converge to 0.");
724                if (info <= A.rows && jobz == 'V')
725                    throw new LapackException("DSYGVD", "Failed to compute an eigenvalue while working on a sub-matrix  " + info + ".");
726                else
727                    throw new LapackException("DSYGVD", "The leading minor of order " + (info - A.rows) + " of B is not positive definite.");
728            }
729        }
730    
731    //END
732    }