001    /*
002     * Copyright (C) 2008-2010 by Holger Arndt
003     *
004     * This file is part of the Universal Java Matrix Package (UJMP).
005     * See the NOTICE file distributed with this work for additional
006     * information regarding copyright ownership and licensing.
007     *
008     * UJMP is free software; you can redistribute it and/or modify
009     * it under the terms of the GNU Lesser General Public License as
010     * published by the Free Software Foundation; either version 2
011     * of the License, or (at your option) any later version.
012     *
013     * UJMP is distributed in the hope that it will be useful,
014     * but WITHOUT ANY WARRANTY; without even the implied warranty of
015     * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
016     * GNU Lesser General Public License for more details.
017     *
018     * You should have received a copy of the GNU Lesser General Public
019     * License along with UJMP; if not, write to the
020     * Free Software Foundation, Inc., 51 Franklin St, Fifth Floor,
021     * Boston, MA  02110-1301  USA
022     */
023    
024    package org.ujmp.core.doublematrix.calculation.general.decomposition;
025    
026    import org.ujmp.core.Matrix;
027    import org.ujmp.core.MatrixFactory;
028    import org.ujmp.core.util.DecompositionOps;
029    import org.ujmp.core.util.MathUtil;
030    import org.ujmp.core.util.UJMPSettings;
031    
032    /**
033     * Singular Value Decomposition.
034     * <P>
035     * For an m-by-n matrix A, the singular value decomposition is an m-by-(m or n)
036     * orthogonal matrix U, a (m or n)-by-n diagonal matrix S, and an n-by-n
037     * orthogonal matrix V so that A = U*S*V'.
038     * <P>
039     * The singular values, sigma[k] = S[k][k], are ordered so that sigma[0] >=
040     * sigma[1] >= ... >= sigma[n-1].
041     * <P>
042     * The singular value decompostion always exists, so the constructor will never
043     * fail. The matrix condition number and the effective numerical rank can be
044     * computed from this decomposition.
045     */
046    
047    public interface SVD<T> {
048    
049            public T[] calc(T source);
050    
051            public static int THRESHOLD = 100;
052    
053            public static final SVD<Matrix> MATRIX = new SVD<Matrix>() {
054    
055                    public final Matrix[] calc(Matrix source) {
056                            if (UJMPSettings.getNumberOfThreads() == 1) {
057                                    if (source.getRowCount() >= THRESHOLD && source.getColumnCount() >= THRESHOLD) {
058                                            return MATRIXLARGESINGLETHREADED.calc(source);
059                                    } else {
060                                            return MATRIXSMALLSINGLETHREADED.calc(source);
061                                    }
062                            } else {
063                                    if (source.getRowCount() >= THRESHOLD && source.getColumnCount() >= THRESHOLD) {
064                                            return MATRIXLARGEMULTITHREADED.calc(source);
065                                    } else {
066                                            return MATRIXSMALLMULTITHREADED.calc(source);
067                                    }
068                            }
069                    }
070            };
071    
072            public static final SVD<Matrix> INSTANCE = MATRIX;
073    
074            public static final SVD<Matrix> UJMP = new SVD<Matrix>() {
075    
076                    public final Matrix[] calc(Matrix source) {
077                            SVDMatrix svd = new SVDMatrix(source);
078                            return new Matrix[] { svd.getU(), svd.getS(), svd.getV() };
079                    }
080            };
081    
082            public static final SVD<Matrix> MATRIXSMALLSINGLETHREADED = UJMP;
083    
084            public static final SVD<Matrix> MATRIXLARGESINGLETHREADED = new SVD<Matrix>() {
085    
086                    public final Matrix[] calc(Matrix source) {
087                            SVD<Matrix> svd = null;
088                            if (UJMPSettings.isUseMTJ()) {
089                                    svd = DecompositionOps.SVD_MTJ;
090                            }
091                            if (svd == null && UJMPSettings.isUseOjalgo()) {
092                                    svd = DecompositionOps.SVD_OJALGO;
093                            }
094                            if (svd == null) {
095                                    svd = UJMP;
096                            }
097                            return svd.calc(source);
098                    }
099            };
100    
101            public static final SVD<Matrix> MATRIXSMALLMULTITHREADED = UJMP;
102    
103            public static final SVD<Matrix> MATRIXLARGEMULTITHREADED = new SVD<Matrix>() {
104    
105                    public final Matrix[] calc(Matrix source) {
106                            SVD<Matrix> svd = null;
107                            if (UJMPSettings.isUseOjalgo()) {
108                                    svd = DecompositionOps.SVD_OJALGO;
109                            }
110                            if (svd == null && UJMPSettings.isUseMTJ()) {
111                                    svd = DecompositionOps.SVD_MTJ;
112                            }
113                            if (svd == null) {
114                                    svd = UJMP;
115                            }
116                            return svd.calc(source);
117                    }
118            };
119    
120            final class SVDMatrix {
121                    private static final double EPSILON = Math.pow(2.0, -52.0);
122    
123                    private static final double TINY = Math.pow(2.0, -966.0);
124    
125                    /**
126                     * Arrays for internal storage of U and V.
127                     * 
128                     * @serial internal storage of U.
129                     * @serial internal storage of V.
130                     */
131                    private final double[][] U, V;
132    
133                    /**
134                     * Array for internal storage of singular values.
135                     * 
136                     * @serial internal storage of singular values.
137                     */
138                    private final double[] s;
139    
140                    /**
141                     * Row and column dimensions.
142                     * 
143                     * @serial row dimension.
144                     * @serial column dimension.
145                     * @serial U column dimension.
146                     */
147                    private final int m, n, ncu;
148    
149                    /**
150                     * Column specification of matrix U
151                     * 
152                     * @serial U column dimension toggle
153                     */
154    
155                    private final boolean thin;
156    
157                    /*
158                     * ------------------------ Old Constructor ------------------------
159                     */
160                    /**
161                     * Construct the singular value decomposition
162                     * 
163                     * @param Arg
164                     *            Rectangular matrix
165                     * @return Structure to access U, S and V.
166                     */
167    
168                    public SVDMatrix(Matrix Arg) {
169                            this(Arg, true, true, true);
170                    }
171    
172                    /*
173                     * ------------------------ Constructor ------------------------
174                     */
175    
176                    /**
177                     * Construct the singular value decomposition
178                     * 
179                     * @param Arg
180                     *            Rectangular matrix
181                     * @param thin
182                     *            If true U is economy sized
183                     * @param wantu
184                     *            If true generate the U matrix
185                     * @param wantv
186                     *            If true generate the V matrix
187                     * @return Structure to access U, S and V.
188                     */
189    
190                    public SVDMatrix(Matrix Arg, boolean thin, boolean wantu, boolean wantv) {
191    
192                            // Derived from LINPACK code.
193                            // Initialize.
194                            final double[][] A = Arg.toDoubleArray();
195                            m = (int) Arg.getRowCount();
196                            n = (int) Arg.getColumnCount();
197                            this.thin = thin;
198    
199                            ncu = thin ? Math.min(m, n) : m;
200                            s = new double[Math.min(m + 1, n)];
201                            U = new double[m][ncu];
202                            V = new double[n][n];
203                            final double[] e = new double[n];
204                            final double[] work = new double[m];
205    
206                            // Reduce A to bidiagonal form, storing the diagonal elements
207                            // in s and the super-diagonal elements in e.
208    
209                            final int nct = Math.min(m - 1, n);
210                            final int nrt = Math.max(0, Math.min(n - 2, m));
211                            final int lu = Math.max(nct, nrt);
212                            for (int k = 0; k < lu; k++) {
213                                    if (k < nct) {
214    
215                                            // Compute the transformation for the k-th column and
216                                            // place the k-th diagonal in s[k].
217                                            // Compute 2-norm of k-th column without under/overflow.
218                                            s[k] = 0;
219                                            for (int i = k; i < m; i++) {
220                                                    s[k] = MathUtil.hypot(s[k], A[i][k]);
221                                            }
222                                            if (s[k] != 0.0) {
223                                                    if (A[k][k] < 0.0) {
224                                                            s[k] = -s[k];
225                                                    }
226                                                    for (int i = k; i < m; i++) {
227                                                            A[i][k] /= s[k];
228                                                    }
229                                                    A[k][k] += 1.0;
230                                            }
231                                            s[k] = -s[k];
232                                    }
233                                    for (int j = k + 1; j < n; j++) {
234                                            if ((k < nct) & (s[k] != 0.0)) {
235    
236                                                    // Apply the transformation.
237    
238                                                    double t = 0;
239                                                    for (int i = k; i < m; i++) {
240                                                            t += A[i][k] * A[i][j];
241                                                    }
242                                                    t = -t / A[k][k];
243                                                    for (int i = k; i < m; i++) {
244                                                            A[i][j] += t * A[i][k];
245                                                    }
246                                            }
247    
248                                            // Place the k-th row of A into e for the
249                                            // subsequent calculation of the row transformation.
250    
251                                            e[j] = A[k][j];
252                                    }
253                                    if (wantu & (k < nct)) {
254    
255                                            // Place the transformation in U for subsequent back
256                                            // multiplication.
257    
258                                            for (int i = k; i < m; i++) {
259                                                    U[i][k] = A[i][k];
260                                            }
261                                    }
262                                    if (k < nrt) {
263    
264                                            // Compute the k-th row transformation and place the
265                                            // k-th super-diagonal in e[k].
266                                            // Compute 2-norm without under/overflow.
267                                            e[k] = 0;
268                                            for (int i = k + 1; i < n; i++) {
269                                                    e[k] = MathUtil.hypot(e[k], e[i]);
270                                            }
271                                            if (e[k] != 0.0) {
272                                                    if (e[k + 1] < 0.0) {
273                                                            e[k] = -e[k];
274                                                    }
275                                                    for (int i = k + 1; i < n; i++) {
276                                                            e[i] /= e[k];
277                                                    }
278                                                    e[k + 1] += 1.0;
279                                            }
280                                            e[k] = -e[k];
281                                            if ((k + 1 < m) & (e[k] != 0.0)) {
282    
283                                                    // Apply the transformation.
284    
285                                                    for (int i = k + 1; i < m; i++) {
286                                                            work[i] = 0.0;
287                                                    }
288                                                    for (int j = k + 1; j < n; j++) {
289                                                            for (int i = k + 1; i < m; i++) {
290                                                                    work[i] += e[j] * A[i][j];
291                                                            }
292                                                    }
293                                                    for (int j = k + 1; j < n; j++) {
294                                                            double t = -e[j] / e[k + 1];
295                                                            for (int i = k + 1; i < m; i++) {
296                                                                    A[i][j] += t * work[i];
297                                                            }
298                                                    }
299                                            }
300                                            if (wantv) {
301    
302                                                    // Place the transformation in V for subsequent
303                                                    // back multiplication.
304    
305                                                    for (int i = k + 1; i < n; i++) {
306                                                            V[i][k] = e[i];
307                                                    }
308                                            }
309                                    }
310                            }
311    
312                            // Set up the final bidiagonal matrix or order p.
313                            int p = Math.min(n, m + 1);
314                            if (nct < n) {
315                                    s[nct] = A[nct][nct];
316                            }
317                            if (m < p) {
318                                    s[p - 1] = 0.0;
319                            }
320                            if (nrt + 1 < p) {
321                                    e[nrt] = A[nrt][p - 1];
322                            }
323                            e[p - 1] = 0.0;
324    
325                            // If required, generate U.
326                            if (wantu) {
327                                    for (int j = nct; j < ncu; j++) {
328                                            for (int i = 0; i < m; i++) {
329                                                    U[i][j] = 0.0;
330                                            }
331                                            U[j][j] = 1.0;
332                                    }
333                                    for (int k = nct - 1; k >= 0; k--) {
334                                            if (s[k] != 0.0) {
335                                                    for (int j = k + 1; j < ncu; j++) {
336                                                            double t = 0;
337                                                            for (int i = k; i < m; i++) {
338                                                                    t += U[i][k] * U[i][j];
339                                                            }
340                                                            t = -t / U[k][k];
341                                                            for (int i = k; i < m; i++) {
342                                                                    U[i][j] += t * U[i][k];
343                                                            }
344                                                    }
345                                                    for (int i = k; i < m; i++) {
346                                                            U[i][k] = -U[i][k];
347                                                    }
348                                                    U[k][k] += 1.0;
349                                                    for (int i = 0; i < k - 1; i++) {
350                                                            U[i][k] = 0.0;
351                                                    }
352                                            } else {
353                                                    for (int i = 0; i < m; i++) {
354                                                            U[i][k] = 0.0;
355                                                    }
356                                                    U[k][k] = 1.0;
357                                            }
358                                    }
359                            }
360    
361                            // If required, generate V.
362                            if (wantv) {
363                                    for (int k = n - 1; k >= 0; k--) {
364                                            if ((k < nrt) & (e[k] != 0.0)) {
365                                                    for (int j = k + 1; j < n; j++) {
366                                                            double t = 0;
367                                                            for (int i = k + 1; i < n; i++) {
368                                                                    t += V[i][k] * V[i][j];
369                                                            }
370                                                            t = -t / V[k + 1][k];
371                                                            for (int i = k + 1; i < n; i++) {
372                                                                    V[i][j] += t * V[i][k];
373                                                            }
374                                                    }
375                                            }
376                                            for (int i = 0; i < n; i++) {
377                                                    V[i][k] = 0.0;
378                                            }
379                                            V[k][k] = 1.0;
380                                    }
381                            }
382    
383                            // Main iteration loop for the singular values.
384    
385                            final int pp = p - 1;
386                            int iter = 0;
387    
388                            while (p > 0) {
389                                    int k, kase;
390    
391                                    // Here is where a test for too many iterations would go.
392    
393                                    // This section of the program inspects for
394                                    // negligible elements in the s and e arrays. On
395                                    // completion the variables kase and k are set as follows.
396    
397                                    // kase = 1 if s(p) and e[k-1] are negligible and k<p
398                                    // kase = 2 if s(k) is negligible and k<p
399                                    // kase = 3 if e[k-1] is negligible, k<p, and
400                                    // s(k), ..., s(p) are not negligible (qr step).
401                                    // kase = 4 if e(p-1) is negligible (convergence).
402    
403                                    for (k = p - 2; k >= -1; k--) {
404                                            if (k == -1) {
405                                                    break;
406                                            }
407                                            if (Math.abs(e[k]) <= TINY + EPSILON * (Math.abs(s[k]) + Math.abs(s[k + 1]))) {
408                                                    e[k] = 0.0;
409                                                    break;
410                                            }
411                                    }
412                                    if (k == p - 2) {
413                                            kase = 4;
414                                    } else {
415                                            int ks;
416                                            for (ks = p - 1; ks >= k; ks--) {
417                                                    if (ks == k) {
418                                                            break;
419                                                    }
420                                                    double t = (ks != p ? Math.abs(e[ks]) : 0.)
421                                                                    + (ks != k + 1 ? Math.abs(e[ks - 1]) : 0.);
422                                                    if (Math.abs(s[ks]) <= TINY + EPSILON * t) {
423                                                            s[ks] = 0.0;
424                                                            break;
425                                                    }
426                                            }
427                                            if (ks == k) {
428                                                    kase = 3;
429                                            } else if (ks == p - 1) {
430                                                    kase = 1;
431                                            } else {
432                                                    kase = 2;
433                                                    k = ks;
434                                            }
435                                    }
436                                    k++;
437    
438                                    // Perform the task indicated by kase.
439    
440                                    switch (kase) {
441    
442                                    // Deflate negligible s(p).
443    
444                                    case 1: {
445                                            double f = e[p - 2];
446                                            e[p - 2] = 0.0;
447                                            for (int j = p - 2; j >= k; j--) {
448                                                    double t = MathUtil.hypot(s[j], f);
449                                                    double cs = s[j] / t;
450                                                    double sn = f / t;
451                                                    s[j] = t;
452                                                    if (j != k) {
453                                                            f = -sn * e[j - 1];
454                                                            e[j - 1] = cs * e[j - 1];
455                                                    }
456                                                    if (wantv) {
457                                                            for (int i = 0; i < n; i++) {
458                                                                    t = cs * V[i][j] + sn * V[i][p - 1];
459                                                                    V[i][p - 1] = -sn * V[i][j] + cs * V[i][p - 1];
460                                                                    V[i][j] = t;
461                                                            }
462                                                    }
463                                            }
464                                    }
465                                            break;
466    
467                                    // Split at negligible s(k).
468    
469                                    case 2: {
470                                            double f = e[k - 1];
471                                            e[k - 1] = 0.0;
472                                            for (int j = k; j < p; j++) {
473                                                    double t = MathUtil.hypot(s[j], f);
474                                                    double cs = s[j] / t;
475                                                    double sn = f / t;
476                                                    s[j] = t;
477                                                    f = -sn * e[j];
478                                                    e[j] = cs * e[j];
479                                                    if (wantu) {
480                                                            for (int i = 0; i < m; i++) {
481                                                                    t = cs * U[i][j] + sn * U[i][k - 1];
482                                                                    U[i][k - 1] = -sn * U[i][j] + cs * U[i][k - 1];
483                                                                    U[i][j] = t;
484                                                            }
485                                                    }
486                                            }
487                                    }
488                                            break;
489    
490                                    // Perform one qr step.
491    
492                                    case 3: {
493    
494                                            // Calculate the shift.
495    
496                                            final double scale = Math.max(Math.max(Math.max(Math.max(Math.abs(s[p - 1]),
497                                                            Math.abs(s[p - 2])), Math.abs(e[p - 2])), Math.abs(s[k])), Math
498                                                            .abs(e[k]));
499                                            final double sp = s[p - 1] / scale;
500                                            final double spm1 = s[p - 2] / scale;
501                                            final double epm1 = e[p - 2] / scale;
502                                            final double sk = s[k] / scale;
503                                            final double ek = e[k] / scale;
504                                            final double b = ((spm1 + sp) * (spm1 - sp) + epm1 * epm1) / 2.0;
505                                            final double c = (sp * epm1) * (sp * epm1);
506                                            double shift = 0.0;
507                                            if ((b != 0.0) | (c != 0.0)) {
508                                                    shift = Math.sqrt(b * b + c);
509                                                    if (b < 0.0) {
510                                                            shift = -shift;
511                                                    }
512                                                    shift = c / (b + shift);
513                                            }
514                                            double f = (sk + sp) * (sk - sp) + shift;
515                                            double g = sk * ek;
516    
517                                            // Chase zeros.
518    
519                                            for (int j = k; j < p - 1; j++) {
520                                                    double t = MathUtil.hypot(f, g);
521                                                    double cs = f / t;
522                                                    double sn = g / t;
523                                                    if (j != k) {
524                                                            e[j - 1] = t;
525                                                    }
526                                                    f = cs * s[j] + sn * e[j];
527                                                    e[j] = cs * e[j] - sn * s[j];
528                                                    g = sn * s[j + 1];
529                                                    s[j + 1] = cs * s[j + 1];
530                                                    if (wantv) {
531                                                            for (int i = 0; i < n; i++) {
532                                                                    t = cs * V[i][j] + sn * V[i][j + 1];
533                                                                    V[i][j + 1] = -sn * V[i][j] + cs * V[i][j + 1];
534                                                                    V[i][j] = t;
535                                                            }
536                                                    }
537                                                    t = MathUtil.hypot(f, g);
538                                                    cs = f / t;
539                                                    sn = g / t;
540                                                    s[j] = t;
541                                                    f = cs * e[j] + sn * s[j + 1];
542                                                    s[j + 1] = -sn * e[j] + cs * s[j + 1];
543                                                    g = sn * e[j + 1];
544                                                    e[j + 1] = cs * e[j + 1];
545                                                    if (wantu && (j < m - 1)) {
546                                                            for (int i = 0; i < m; i++) {
547                                                                    t = cs * U[i][j] + sn * U[i][j + 1];
548                                                                    U[i][j + 1] = -sn * U[i][j] + cs * U[i][j + 1];
549                                                                    U[i][j] = t;
550                                                            }
551                                                    }
552                                            }
553                                            e[p - 2] = f;
554                                            iter++;
555                                    }
556                                            break;
557    
558                                    // Convergence.
559    
560                                    case 4: {
561    
562                                            // Make the singular values positive.
563    
564                                            if (s[k] <= 0.0) {
565                                                    s[k] = (s[k] < 0.0 ? -s[k] : 0.0);
566                                                    if (wantv) {
567                                                            for (int i = 0; i < n; i++) {
568                                                                    V[i][k] = -V[i][k];
569                                                            }
570                                                    }
571                                            }
572    
573                                            // Order the singular values.
574    
575                                            while (k < pp) {
576                                                    if (s[k] >= s[k + 1]) {
577                                                            break;
578                                                    }
579                                                    double t = s[k];
580                                                    s[k] = s[k + 1];
581                                                    s[k + 1] = t;
582                                                    if (wantv && (k < n - 1)) {
583                                                            for (int i = 0; i < n; i++) {
584                                                                    t = V[i][k + 1];
585                                                                    V[i][k + 1] = V[i][k];
586                                                                    V[i][k] = t;
587                                                            }
588                                                    }
589                                                    if (wantu && (k < m - 1)) {
590                                                            for (int i = 0; i < m; i++) {
591                                                                    t = U[i][k + 1];
592                                                                    U[i][k + 1] = U[i][k];
593                                                                    U[i][k] = t;
594                                                            }
595                                                    }
596                                                    k++;
597                                            }
598                                            iter = 0;
599                                            p--;
600                                    }
601                                            break;
602                                    }
603                            }
604                    }
605    
606                    /*
607                     * ------------------------ Public Methods ------------------------
608                     */
609    
610                    /**
611                     * Return the left singular vectors
612                     * 
613                     * @return U
614                     */
615    
616                    public final Matrix getU() {
617                            final double[][] x = new double[m][m >= n ? (thin ? Math.min(m + 1, n) : ncu) : ncu];
618    
619                            for (int r = 0; r < m; r++) {
620                                    for (int c = x[0].length; --c >= 0;) {
621                                            x[r][c] = U[r][c];
622                                    }
623                            }
624    
625                            return MatrixFactory.linkToArray(x);
626                    }
627    
628                    /**
629                     * Return the right singular vectors
630                     * 
631                     * @return V
632                     */
633    
634                    public final Matrix getV() {
635                            return V == null ? null : MatrixFactory.linkToArray(V);
636                    }
637    
638                    /**
639                     * Return the one-dimensional array of singular values
640                     * 
641                     * @return diagonal of S.
642                     */
643    
644                    public final double[] getSingularValues() {
645                            return s;
646                    }
647    
648                    /**
649                     * Return the diagonal matrix of singular values
650                     * 
651                     * @return S
652                     */
653    
654                    public final Matrix getS() {
655                            final double[][] X = new double[m >= n ? (thin ? n : ncu) : ncu][n];
656                            for (int i = Math.min(m, n); --i >= 0;)
657                                    X[i][i] = s[i];
658                            return MatrixFactory.linkToArray(X);
659                    }
660    
661                    /**
662                     * Return the diagonal matrix of the reciprocals of the singular values
663                     * 
664                     * @return S+
665                     */
666    
667                    public final Matrix getreciprocalS() {
668                            final double[][] X = new double[n][m >= n ? (thin ? n : ncu) : ncu];
669                            for (int i = Math.min(m, n) - 1; i >= 0; i--)
670                                    X[i][i] = s[i] == 0.0 ? 0.0 : 1.0 / s[i];
671                            return MatrixFactory.linkToArray(X);
672                    }
673    
674                    /**
675                     * Return the Moore-Penrose (generalized) inverse Slightly modified
676                     * version of Kim van der Linde's code
677                     * 
678                     * @param omit
679                     *            if true tolerance based omitting of negligible singular
680                     *            values
681                     * @return A+
682                     */
683    
684                    public final Matrix inverse(boolean omit) {
685                            final double[][] inverse = new double[n][m];
686                            if (rank() > 0) {
687                                    final double[] reciprocalS = new double[s.length];
688                                    if (omit) {
689                                            double tol = Math.max(m, n) * s[0] * EPSILON;
690                                            for (int i = s.length - 1; i >= 0; i--)
691                                                    reciprocalS[i] = Math.abs(s[i]) < tol ? 0.0 : 1.0 / s[i];
692                                    } else
693                                            for (int i = s.length - 1; i >= 0; i--)
694                                                    reciprocalS[i] = s[i] == 0.0 ? 0.0 : 1.0 / s[i];
695                                    int min = Math.min(n, ncu);
696                                    for (int i = n - 1; i >= 0; i--)
697                                            for (int j = m - 1; j >= 0; j--)
698                                                    for (int k = min - 1; k >= 0; k--)
699                                                            inverse[i][j] += V[i][k] * reciprocalS[k] * U[j][k];
700                            }
701                            return MatrixFactory.linkToArray(inverse);
702                    }
703    
704                    /**
705                     * Two norm
706                     * 
707                     * @return max(S)
708                     */
709    
710                    public final double norm2() {
711                            return s[0];
712                    }
713    
714                    /**
715                     * Two norm condition number
716                     * 
717                     * @return max(S)/min(S)
718                     */
719    
720                    public final double cond() {
721                            return s[0] / s[Math.min(m, n) - 1];
722                    }
723    
724                    /**
725                     * Effective numerical matrix rank
726                     * 
727                     * @return Number of nonnegligible singular values.
728                     */
729    
730                    public final int rank() {
731                            final double tol = Math.max(m, n) * s[0] * EPSILON;
732                            int r = 0;
733                            for (int i = 0; i < s.length; i++) {
734                                    if (s[i] > tol) {
735                                            r++;
736                                    }
737                            }
738                            return r;
739                    }
740    
741            }
742    }