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