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.doublematrix.DenseDoubleMatrix2D;
028    import org.ujmp.core.doublematrix.calculation.AbstractDoubleCalculation;
029    import org.ujmp.core.doublematrix.impl.ArrayDenseDoubleMatrix2D;
030    import org.ujmp.core.exceptions.MatrixException;
031    import org.ujmp.core.interfaces.HasRowMajorDoubleArray2D;
032    import org.ujmp.core.util.UJMPSettings;
033    
034    /**
035     * <p>
036     * This class implements some matrix operations that I need for other things I'm
037     * doing. Why not use existing packages? Well, I've tried, but they seem to not
038     * have simple linear algebra things like inverting singular or non-square
039     * inverses. Maybe I've just been looking in the wrong places.
040     * </p>
041     * <p>
042     * Anyway I wrote this algorithm in 1989 using Ada and Fortran77. In 1992 I
043     * converted it to C++, and in 2003 I made it a C++ template to be used for
044     * absolutely anything that you'd want.
045     * </p>
046     * <p>
047     * I attempted to convert this to a Java template, but I need to make objects of
048     * the parameter type, and without doing some inspection or copy tricks, I can't
049     * get there from here. I only use this for doubles anyway, so here's the
050     * version that I use. The C++ version is still available.
051     * <p>
052     * </p>
053     * The matrix inversion was in from the start, as the only really useful part of
054     * the Class. I added the bandwidth reduction routine in 1991 - I was stuck in
055     * SOS (USAF school) at the time and was thinking about optimizing the bandwidth
056     * of a matrix made from a finite-element grid by renumbering the nodes of the
057     * grid. </p>
058     * <p>
059     * Changes by Holger Arndt: The original code has been adapted for the Universal
060     * Java Matrix Package. Methods for different matrix implementations have been
061     * added.
062     * </p>
063     * 
064     * @author Rand Huso
065     * @author Holger Arndt
066     */
067    public class Ginv extends AbstractDoubleCalculation {
068            private static final long serialVersionUID = 1087531579133023922L;
069    
070            public Ginv(Matrix source) {
071                    super(source);
072            }
073    
074            /**
075             * This routine performs the matrix multiplication. The final matrix size is
076             * taken from the rows of the left matrix and the columns of the right
077             * matrix. The timesInner is the minimum of the left columns and the right
078             * rows.
079             * 
080             * @param matrix1
081             *            the first matrix
082             * @param matrix2
083             *            the second matrix
084             * @param timesInner
085             *            number of rows/columns to process
086             * @return product of the two matrices
087             */
088            public static DenseDoubleMatrix2D times(DenseDoubleMatrix2D matrix1,
089                            DenseDoubleMatrix2D matrix2, long timesInner) {
090                    long timesRows = matrix1.getRowCount();
091                    long timesCols = matrix2.getColumnCount();
092                    DenseDoubleMatrix2D response = DenseDoubleMatrix2D.factory.zeros(timesRows, timesCols);
093                    for (long row = 0; row < timesRows; row++) {
094                            for (long col = 0; col < timesCols; col++) {
095                                    for (long inner = 0; inner < timesInner; inner++) {
096                                            response.setDouble(matrix1.getAsDouble(row, inner)
097                                                            * matrix2.getDouble(inner, col) + response.getDouble(row, col), row,
098                                                            col);
099                                    }
100                            }
101                    }
102                    return response;
103            }
104    
105            /**
106             * This routine performs the matrix multiplication. The final matrix size is
107             * taken from the rows of the left matrix and the columns of the right
108             * matrix. The timesInner is the minimum of the left columns and the right
109             * rows.
110             * 
111             * @param matrix1
112             *            the first matrix
113             * @param matrix2
114             *            the second matrix
115             * @param timesInner
116             *            number of rows/columns to process
117             * @return product of the two matrices
118             */
119            public static double[][] times(double[][] matrix1, double[][] matrix2, int timesInner) {
120                    int timesRows = matrix1.length;
121                    int timesCols = matrix2[0].length;
122                    double[][] response = new double[timesRows][timesCols];
123                    for (int row = 0; row < timesRows; row++) {
124                            for (int col = 0; col < timesCols; col++) {
125                                    for (int inner = 0; inner < timesInner; inner++) {
126                                            response[row][col] = matrix1[row][inner] * matrix2[inner][col]
127                                                            + response[row][col];
128                                    }
129                            }
130                    }
131                    return response;
132            }
133    
134            /**
135             * Swap components in the two columns.
136             * 
137             * @param matrix
138             *            the matrix to modify
139             * @param col1
140             *            the first row
141             * @param col2
142             *            the second row
143             */
144            public static void swapCols(Matrix matrix, long col1, long col2) {
145                    double temp = 0;
146                    long rows = matrix.getRowCount();
147                    for (long row = 0; row < rows; row++) {
148                            temp = matrix.getAsDouble(row, col1);
149                            matrix.setAsDouble(matrix.getAsDouble(row, col2), row, col1);
150                            matrix.setAsDouble(temp, row, col2);
151                    }
152            }
153    
154            /**
155             * Swap components in the two columns.
156             * 
157             * @param matrix
158             *            the matrix to modify
159             * @param col1
160             *            the first row
161             * @param col2
162             *            the second row
163             */
164            public static void swapCols(DenseDoubleMatrix2D matrix, long col1, long col2) {
165                    double temp = 0;
166                    long rows = matrix.getRowCount();
167                    for (long row = 0; row < rows; row++) {
168                            temp = matrix.getDouble(row, col1);
169                            matrix.setDouble(matrix.getDouble(row, col2), row, col1);
170                            matrix.setDouble(temp, row, col2);
171                    }
172            }
173    
174            /**
175             * Swap components in the two columns.
176             * 
177             * @param matrix
178             *            the matrix to modify
179             * @param col1
180             *            the first row
181             * @param col2
182             *            the second row
183             */
184            public static void swapCols(double[][] matrix, int col1, int col2) {
185                    double temp = 0;
186                    int rows = matrix.length;
187                    double[] r = null;
188                    for (int row = 0; row < rows; row++) {
189                            r = matrix[row];
190                            temp = r[col1];
191                            r[col1] = r[col2];
192                            r[col2] = temp;
193                    }
194            }
195    
196            /**
197             * Swap components in the two rows.
198             * 
199             * @param matrix
200             *            the matrix to modify
201             * @param row1
202             *            the first row
203             * @param row2
204             *            the second row
205             */
206            public static void swapRows(Matrix matrix, long row1, long row2) {
207                    double temp = 0;
208                    long cols = matrix.getColumnCount();
209                    for (long col = 0; col < cols; col++) {
210                            temp = matrix.getAsDouble(row1, col);
211                            matrix.setAsDouble(matrix.getAsDouble(row2, col), row1, col);
212                            matrix.setAsDouble(temp, row2, col);
213                    }
214            }
215    
216            /**
217             * Swap components in the two rows.
218             * 
219             * @param matrix
220             *            the matrix to modify
221             * @param row1
222             *            the first row
223             * @param row2
224             *            the second row
225             */
226            public static void swapRows(DenseDoubleMatrix2D matrix, long row1, long row2) {
227                    double temp = 0;
228                    long cols = matrix.getColumnCount();
229                    for (long col = 0; col < cols; col++) {
230                            temp = matrix.getDouble(row1, col);
231                            matrix.setDouble(matrix.getDouble(row2, col), row1, col);
232                            matrix.setDouble(temp, row2, col);
233                    }
234            }
235    
236            /**
237             * Swap components in the two rows.
238             * 
239             * @param matrix
240             *            the matrix to modify
241             * @param row1
242             *            the first row
243             * @param row2
244             *            the second row
245             */
246            public static void swapRows(double[][] matrix, int row1, int row2) {
247                    double[] temp = matrix[row1];
248                    matrix[row1] = matrix[row2];
249                    matrix[row2] = temp;
250            }
251    
252            /**
253             * <p>
254             * Matrix inversion is the reason this Class exists - this method creates a
255             * generalized matrix inverse of the current matrix. The result is returned
256             * in a new matrix.
257             * </p>
258             * <p>
259             * Matrices may be square, non-square, or singular. The operations will be
260             * identical. If the matrix has a single possible inverse, there will be no
261             * arbitrariness to the solution. The method here was suggested to me by
262             * John Jones Jr. at AFIT in the 1980s, and this algorithm is an original
263             * creation of mine to implement his method.
264             * </p>
265             * <p>
266             * A matrix inverse has some properties described here. Let
267             * <code><b>A</b></code> be the original matrix. Let the inverse, as
268             * calculated here be <code><b>A12</b></code> (an inverse with properties 1
269             * and 2 - being left side inverse and right side inverse). An inverse times
270             * the original matrix should yield an identity matrix.
271             * <code><b>A x = b</b></code> is the usual equation where
272             * <code><b>A</b></code> is the matrix, <code><b>x</b></code> is a vector of
273             * the unknowns and <code><b>b</b></code> is a vector of the constant
274             * values:
275             * </p>
276             * <p>
277             * Given these equations:
278             * <code><b>C x + D y + E z = b1 F x + G y + H z = b1</b></code>
279             * </p>
280             * <p>
281             * (The usual programs available don't handle more unknowns than equations,
282             * and will stop at this point)
283             * </p>
284             * <p>
285             * The <code><b>A</b></code> matrix is:
286             * <code><b>| C D E | | F G H |</b></code>
287             * </p>
288             * <p>
289             * The <code><b>x</b></code> vector is:
290             * <code><b>| x | | y | | z |</b></code>
291             * </p>
292             * <p>
293             * The <code><b>b</b></code> vector is: <code><b>| b1 | | b2 |</b></code>
294             * </p>
295             * <p>
296             * <code><b>A * x = b</b></code>
297             * </p>
298             * <p>
299             * The generalized inverse <code><b>A12</b></code> in this case will be of
300             * size (3,2): <code><b>| J K | | L M | | N P |</b></code>
301             * </p>
302             * <p>
303             * The left-hand inverse is defined such that the product of the
304             * (generalized) inverse times the original matrix times the (generalized)
305             * inverse will yield the (generalized) inverse again:
306             * <code><b>A12 * (A * A12) = A12</b></code>
307             * </p>
308             * <p>
309             * The right-hand inverse is defined similarly:
310             * <code><b>(A * A12) * A = A</b></code>
311             * </p>
312             * <p>
313             * If a matrix (<code><b>A12</b></code>) meets these criteria, it's
314             * considered to be a generalized inverse of the <code><b>A</b></code>
315             * matrix, even though it may not be square, or the product of
316             * <code><b>A * A12</b></code> or <code><b>A12 * A</b></code> may not be the
317             * identity matrix! (Won't be if the input <code><b>A</b></code> matrix is
318             * non-square or singular)
319             * </p>
320             * <p>
321             * In the case of <code><b>(A * A12)</b></code> being the identity matrix,
322             * the product of <code><b>(A12 *
323             * A)</b></code> will also be the identity matrix, and the solution will be
324             * unique: <code><b>A12</b></code> will be the exact and only solution to
325             * the equation.
326             * </p>
327             * <p>
328             * Refer to {@link http://mjollnir.com/matrix/} for the best description.
329             * </p>
330             * 
331             * @param matrix
332             *            matrix to invert
333             * @return a generalized matrix inverse (possibly not unique)
334             */
335            public static DenseDoubleMatrix2D inverse(Matrix matrix) {
336                    double epsilon = UJMPSettings.getTolerance();
337                    long rows = matrix.getRowCount();
338                    long cols = matrix.getColumnCount();
339                    DenseDoubleMatrix2D s = DenseDoubleMatrix2D.factory.zeros(cols, cols);
340                    s.eye(Ret.ORIG);
341                    DenseDoubleMatrix2D t = DenseDoubleMatrix2D.factory.zeros(rows, rows);
342                    t.eye(Ret.ORIG);
343                    long maxDiag = Math.min(rows, cols);
344    
345                    int diag = 0;
346                    for (; diag < maxDiag; diag++) {
347    
348                            // get the largest value for the pivot
349                            swapPivot(matrix, diag, s, t);
350    
351                            if (matrix.getAsDouble(diag, diag) == 0.0) {
352                                    break;
353                            }
354    
355                            // divide through to make pivot identity
356                            double divisor = matrix.getAsDouble(diag, diag);
357                            if (Math.abs(divisor) < epsilon) {
358                                    matrix.setAsDouble(0.0, diag, diag);
359                                    break;
360                            }
361    
362                            divideRowBy(matrix, diag, diag, divisor);
363                            divideRowBy(t, diag, 0, divisor);
364                            matrix.setAsDouble(1.0, diag, diag);
365    
366                            // remove values down remaining rows
367                            for (long row = diag + 1; row < rows; row++) {
368                                    double factor = matrix.getAsDouble(row, diag);
369                                    if (factor != 0.0) {
370                                            addRowTimes(matrix, diag, diag, row, factor);
371                                            addRowTimes(t, diag, 0, row, factor);
372                                            matrix.setAsDouble(0.0, row, diag);
373                                    }
374                            }
375    
376                            // remove values across remaining cols - some optimization could
377                            // be done here because the changes to the original matrix at this
378                            // point only touch the current diag column
379                            for (long col = diag + 1; col < cols; col++) {
380                                    double factor = matrix.getAsDouble(diag, col);
381                                    if (factor != 0.0) {
382                                            addColTimes(matrix, diag, diag, col, factor);
383                                            addColTimes(s, diag, 0, col, factor);
384                                            matrix.setAsDouble(0.0, diag, col);
385                                    }
386                            }
387                    }
388    
389                    return times(s, t, diag);
390            }
391    
392            /**
393             * Same as {@link inverse(Matrix)} but optimized for 2D dense double
394             * matrices
395             * 
396             * @param matrix
397             *            the matrix to invert
398             * @return generalized matrix inverse
399             */
400            public static DenseDoubleMatrix2D inverse(DenseDoubleMatrix2D matrix) {
401                    double epsilon = UJMPSettings.getTolerance();
402                    long rows = matrix.getRowCount();
403                    long cols = matrix.getColumnCount();
404                    DenseDoubleMatrix2D s = DenseDoubleMatrix2D.factory.zeros(cols, cols);
405                    s.eye(Ret.ORIG);
406                    DenseDoubleMatrix2D t = DenseDoubleMatrix2D.factory.zeros(rows, rows);
407                    t.eye(Ret.ORIG);
408                    long maxDiag = Math.min(rows, cols);
409    
410                    int diag = 0;
411                    for (; diag < maxDiag; diag++) {
412    
413                            // get the largest value for the pivot
414                            swapPivot(matrix, diag, s, t);
415    
416                            if (matrix.getAsDouble(diag, diag) == 0.0) {
417                                    break;
418                            }
419    
420                            // divide through to make pivot identity
421                            double divisor = matrix.getAsDouble(diag, diag);
422                            if (Math.abs(divisor) < epsilon) {
423                                    matrix.setDouble(0.0, diag, diag);
424                                    break;
425                            }
426    
427                            divideRowBy(matrix, diag, diag, divisor);
428                            divideRowBy(t, diag, 0, divisor);
429                            matrix.setDouble(1.0, diag, diag);
430    
431                            // remove values down remaining rows
432                            for (long row = diag + 1; row < rows; row++) {
433                                    double factor = matrix.getDouble(row, diag);
434                                    if (factor != 0.0) {
435                                            addRowTimes(matrix, diag, diag, row, factor);
436                                            addRowTimes(t, diag, 0, row, factor);
437                                            matrix.setDouble(0.0, row, diag);
438                                    }
439                            }
440    
441                            // remove values across remaining cols - some optimization could
442                            // be done here because the changes to the original matrix at this
443                            // point only touch the current diag column
444                            for (long col = diag + 1; col < cols; col++) {
445                                    double factor = matrix.getDouble(diag, col);
446                                    if (factor != 0.0) {
447                                            addColTimes(matrix, diag, diag, col, factor);
448                                            addColTimes(s, diag, 0, col, factor);
449                                            matrix.setDouble(0.0, diag, col);
450                                    }
451                            }
452                    }
453    
454                    return times(s, t, diag);
455            }
456    
457            /**
458             * Same as {@link inverse(Matrix)} but optimized for 2D double arrays
459             * 
460             * @param matrix
461             *            the matrix to invert
462             * @return generalized matrix inverse
463             */
464            public static DenseDoubleMatrix2D inverse(final double[][] matrix) {
465                    double epsilon = UJMPSettings.getTolerance();
466                    int rows = matrix.length;
467                    int cols = matrix[0].length;
468                    double[][] s = new double[cols][cols];
469                    for (int c = 0; c < cols; c++) {
470                            s[c][c] = 1.0;
471                    }
472                    final double[][] t = new double[rows][rows];
473                    for (int r = 0; r < rows; r++) {
474                            t[r][r] = 1.0;
475                    }
476                    int maxDiag = Math.min(rows, cols);
477    
478                    int diag = 0;
479                    for (; diag < maxDiag; diag++) {
480    
481                            // get the largest value for the pivot
482                            swapPivot(matrix, diag, s, t);
483    
484                            if (matrix[diag][diag] == 0.0) {
485                                    break;
486                            }
487    
488                            // divide through to make pivot identity
489                            double divisor = matrix[diag][diag];
490                            if (Math.abs(divisor) < epsilon) {
491                                    matrix[diag][diag] = 0.0;
492                                    break;
493                            }
494    
495                            divideRowBy(matrix, diag, diag, divisor);
496                            divideRowBy(t, diag, 0, divisor);
497                            matrix[diag][diag] = 1.0;
498    
499                            // remove values down remaining rows
500                            for (int row = diag + 1; row < rows; row++) {
501                                    double factor = matrix[row][diag];
502                                    if (factor != 0.0) {
503                                            addRowTimes(matrix, diag, diag, row, factor);
504                                            addRowTimes(t, diag, 0, row, factor);
505                                            matrix[row][diag] = 0.0;
506                                    }
507                            }
508    
509                            // remove values across remaining cols - some optimization could
510                            // be done here because the changes to the original matrix at this
511                            // point only touch the current diag column
512                            for (int col = diag + 1; col < cols; col++) {
513                                    double factor = matrix[diag][col];
514                                    if (factor != 0.0) {
515                                            addColTimes(matrix, diag, diag, col, factor);
516                                            addColTimes(s, diag, 0, col, factor);
517                                            matrix[diag][col] = 0.0;
518                                    }
519                            }
520                    }
521    
522                    double[][] result = times(s, t, diag);
523                    return new ArrayDenseDoubleMatrix2D(result);
524            }
525    
526            /**
527             * Add a factor times one column to another column
528             * 
529             * @param matrix
530             *            the matrix to modify
531             * @param diag
532             *            coordinate on the diagonal
533             * @param fromRow
534             *            first row to process
535             * @param col
536             *            column to process
537             * @param factor
538             *            factor to multiply
539             */
540            public static void addColTimes(Matrix matrix, long diag, long fromRow, long col, double factor) {
541                    long rows = matrix.getRowCount();
542                    for (long row = fromRow; row < rows; row++) {
543                            matrix.setAsDouble(matrix.getAsDouble(row, col) - factor
544                                            * matrix.getAsDouble(row, diag), row, col);
545                    }
546            }
547    
548            /**
549             * Add a factor times one column to another column
550             * 
551             * @param matrix
552             *            the matrix to modify
553             * @param diag
554             *            coordinate on the diagonal
555             * @param fromRow
556             *            first row to process
557             * @param col
558             *            column to process
559             * @param factor
560             *            factor to multiply
561             */
562            public static void addColTimes(double[][] matrix, int diag, int fromRow, int col, double factor) {
563                    int rows = matrix.length;
564                    double[] r = null;
565                    for (int row = fromRow; row < rows; row++) {
566                            r = matrix[row];
567                            r[col] -= factor * r[diag];
568                    }
569            }
570    
571            /**
572             * Add a factor times one column to another column
573             * 
574             * @param matrix
575             *            the matrix to modify
576             * @param diag
577             *            coordinate on the diagonal
578             * @param fromRow
579             *            first row to process
580             * @param col
581             *            column to process
582             * @param factor
583             *            factor to multiply
584             */
585            public static void addColTimes(DenseDoubleMatrix2D matrix, long diag, long fromRow, long col,
586                            double factor) {
587                    long rows = matrix.getRowCount();
588                    for (long row = fromRow; row < rows; row++) {
589                            matrix.setDouble(matrix.getDouble(row, col) - factor * matrix.getDouble(row, diag),
590                                            row, col);
591                    }
592            }
593    
594            /**
595             * Add a factor times one row to another row
596             * 
597             * @param matrix
598             *            the matrix to modify
599             * @param diag
600             *            coordinate on the diagonal
601             * @param fromCol
602             *            first column to process
603             * @param row
604             *            row to process
605             * @param factor
606             *            factor to multiply
607             */
608            public static void addRowTimes(DenseDoubleMatrix2D matrix, long diag, long fromCol, long row,
609                            double factor) {
610                    long cols = matrix.getColumnCount();
611                    for (long col = fromCol; col < cols; col++) {
612                            matrix.setDouble(matrix.getDouble(row, col) - factor * matrix.getDouble(diag, col),
613                                            row, col);
614                    }
615            }
616    
617            /**
618             * Add a factor times one row to another row
619             * 
620             * @param matrix
621             *            the matrix to modify
622             * @param diag
623             *            coordinate on the diagonal
624             * @param fromCol
625             *            first column to process
626             * @param row
627             *            row to process
628             * @param factor
629             *            factor to multiply
630             */
631            public static void addRowTimes(double[][] matrix, int diag, int fromCol, int row, double factor) {
632                    int cols = matrix[0].length;
633                    double[] d = matrix[diag];
634                    double[] r = matrix[row];
635                    for (int col = fromCol; col < cols; col++) {
636                            r[col] -= factor * d[col];
637                    }
638            }
639    
640            /**
641             * Add a factor times one row to another row
642             * 
643             * @param matrix
644             *            the matrix to modify
645             * @param diag
646             *            coordinate on the diagonal
647             * @param fromCol
648             *            first column to process
649             * @param row
650             *            row to process
651             * @param factor
652             *            factor to multiply
653             */
654            public static void addRowTimes(Matrix matrix, long diag, long fromCol, long row, double factor) {
655                    long cols = matrix.getColumnCount();
656                    for (long col = fromCol; col < cols; col++) {
657                            matrix.setAsDouble(matrix.getAsDouble(row, col) - factor
658                                            * matrix.getAsDouble(diag, col), row, col);
659                    }
660            }
661    
662            /**
663             * Divide the row from this column position by this value
664             * 
665             * @param matrix
666             *            the matrix to modify
667             * @param aRow
668             *            the row to process
669             * @param fromCol
670             *            starting from column
671             * @param value
672             *            the value to divide
673             */
674            public static void divideRowBy(Matrix matrix, long aRow, long fromCol, double value) {
675                    long cols = matrix.getColumnCount();
676                    for (long col = fromCol; col < cols; col++) {
677                            matrix.setAsDouble(matrix.getAsDouble(aRow, col) / value, aRow, col);
678                    }
679            }
680    
681            /**
682             * Divide the row from this column position by this value
683             * 
684             * @param matrix
685             *            the matrix to modify
686             * @param aRow
687             *            the row to process
688             * @param fromCol
689             *            starting from column
690             * @param value
691             *            the value to divide
692             */
693            public static void divideRowBy(DenseDoubleMatrix2D matrix, long aRow, long fromCol, double value) {
694                    long cols = matrix.getColumnCount();
695                    for (long col = fromCol; col < cols; col++) {
696                            matrix.setDouble(matrix.getDouble(aRow, col) / value, aRow, col);
697                    }
698            }
699    
700            /**
701             * Divide the row from this column position by this value
702             * 
703             * @param matrix
704             *            the matrix to modify
705             * @param aRow
706             *            the row to process
707             * @param fromCol
708             *            starting from column
709             * @param value
710             *            the value to divide
711             */
712            public static void divideRowBy(double[][] matrix, int aRow, int fromCol, double value) {
713                    int cols = matrix[0].length;
714                    double[] r = matrix[aRow];
715                    for (int col = fromCol; col < cols; col++) {
716                            r[col] /= value;
717                    }
718            }
719    
720            /**
721             * Swap the matrices so that the largest value is on the pivot
722             * 
723             * @param source
724             *            the matrix to modify
725             * @param diag
726             *            the position on the diagonal
727             * @param s
728             *            the matrix s
729             * @param t
730             *            the matrix t
731             */
732            public static void swapPivot(Matrix source, long diag, Matrix s, Matrix t) {
733                    // get swap row and col
734                    long swapRow = diag;
735                    long swapCol = diag;
736                    double maxValue = Math.abs(source.getAsDouble(diag, diag));
737                    long rows = source.getRowCount();
738                    long cols = source.getColumnCount();
739                    double abs = 0;
740    
741                    for (long row = diag; row < rows; row++) {
742                            for (long col = diag; col < cols; col++) {
743                                    abs = Math.abs(source.getAsDouble(row, col));
744                                    if (abs > maxValue) {
745                                            maxValue = abs;
746                                            swapRow = row;
747                                            swapCol = col;
748                                    }
749                            }
750                    }
751    
752                    // swap rows and columns
753                    if (swapRow != diag) {
754                            swapRows(source, swapRow, diag);
755                            swapRows(t, swapRow, diag);
756                    }
757                    if (swapCol != diag) {
758                            swapCols(source, swapCol, diag);
759                            swapCols(s, swapCol, diag);
760                    }
761            }
762    
763            /**
764             * Swap the matrices so that the largest value is on the pivot
765             * 
766             * @param source
767             *            the matrix to modify
768             * @param diag
769             *            the position on the diagonal
770             * @param s
771             *            the matrix s
772             * @param t
773             *            the matrix t
774             */
775            public static void swapPivot(DenseDoubleMatrix2D source, long diag, DenseDoubleMatrix2D s,
776                            DenseDoubleMatrix2D t) {
777                    // get swap row and col
778                    long swapRow = diag;
779                    long swapCol = diag;
780                    double maxValue = Math.abs(source.getDouble(diag, diag));
781                    long rows = source.getRowCount();
782                    long cols = source.getColumnCount();
783                    double abs = 0;
784                    for (long row = diag; row < rows; row++) {
785                            for (long col = diag; col < cols; col++) {
786                                    abs = Math.abs(source.getDouble(row, col));
787                                    if (abs > maxValue) {
788                                            maxValue = abs;
789                                            swapRow = row;
790                                            swapCol = col;
791                                    }
792                            }
793                    }
794    
795                    // swap rows and columns
796                    if (swapRow != diag) {
797                            swapRows(source, swapRow, diag);
798                            swapRows(t, swapRow, diag);
799                    }
800                    if (swapCol != diag) {
801                            swapCols(source, swapCol, diag);
802                            swapCols(s, swapCol, diag);
803                    }
804            }
805    
806            /**
807             * Swap the matrices so that the largest value is on the pivot
808             * 
809             * @param source
810             *            the matrix to modify
811             * @param diag
812             *            the position on the diagonal
813             * @param s
814             *            the matrix s
815             * @param t
816             *            the matrix t
817             */
818            public static void swapPivot(double[][] source, int diag, double[][] s, double[][] t) {
819                    // get swap row and col
820                    int swapRow = diag;
821                    int swapCol = diag;
822                    double maxValue = Math.abs(source[diag][diag]);
823                    int rows = source.length;
824                    int cols = source[0].length;
825                    double abs = 0;
826                    double[] r = null;
827                    for (int row = diag; row < rows; row++) {
828                            r = source[row];
829                            for (int col = diag; col < cols; col++) {
830                                    abs = Math.abs(r[col]);
831                                    if (abs > maxValue) {
832                                            maxValue = abs;
833                                            swapRow = row;
834                                            swapCol = col;
835                                    }
836                            }
837                    }
838    
839                    // swap rows and columns
840                    if (swapRow != diag) {
841                            swapRows(source, swapRow, diag);
842                            swapRows(t, swapRow, diag);
843                    }
844                    if (swapCol != diag) {
845                            swapCols(source, swapCol, diag);
846                            swapCols(s, swapCol, diag);
847                    }
848            }
849    
850            /**
851             * Check to see if a non-zero and a zero value in the rows leading up to
852             * this column can be swapped. This is part of the bandwidth reduction
853             * algorithm.
854             * 
855             * @param matrix
856             *            the matrix to check
857             * @param row1
858             *            the first row
859             * @param row2
860             *            the second row
861             * @param col1
862             *            the column
863             * @return true if the rows can be swapped
864             */
865            public static boolean canSwapRows(Matrix matrix, int row1, int row2, int col1) {
866                    boolean response = true;
867                    for (int col = 0; col < col1; ++col) {
868                            if (0 == matrix.getAsDouble(row1, col)) {
869                                    if (0 != matrix.getAsDouble(row2, col)) {
870                                            response = false;
871                                            break;
872                                    }
873                            }
874                    }
875                    return response;
876            }
877    
878            /**
879             * Check to see if columns can be swapped - part of the bandwidth reduction
880             * algorithm.
881             * 
882             * @param matrix
883             *            the matrix to check
884             * @param col1
885             *            the first column
886             * @param col2
887             *            the second column
888             * @param row1
889             *            the row
890             * @return true if the columns can be swapped
891             */
892            public static boolean canSwapCols(Matrix matrix, int col1, int col2, int row1) {
893                    boolean response = true;
894                    for (int row = row1 + 1; row < matrix.getRowCount(); ++row) {
895                            if (0 == matrix.getAsDouble(row, col1)) {
896                                    if (0 != matrix.getAsDouble(row, col2)) {
897                                            response = false;
898                                            break;
899                                    }
900                            }
901                    }
902                    return response;
903            }
904    
905            public static Matrix reduce(Matrix source, Matrix response) {
906                    if (source.getRowCount() == source.getColumnCount()) {
907                            // pass one (descending the diagonal):
908                            for (int col = 0; col < source.getColumnCount() - 1; ++col) {
909                                    for (int rowData = (int) source.getRowCount() - 1; rowData > col; --rowData) {
910                                            if (0 != source.getAsDouble(rowData, col)) {
911                                                    for (int rowEmpty = rowData - 1; rowEmpty > col; --rowEmpty) {
912                                                            if (0 == source.getAsDouble(rowEmpty, col)) {
913                                                                    if (Ginv.canSwapRows(source, rowData, rowEmpty, col)) {
914                                                                            Ginv.swapRows(source, rowData, rowEmpty);
915                                                                            Ginv.swapCols(source, rowData, rowEmpty);
916                                                                            Ginv.swapRows(response, rowData, rowEmpty);
917                                                                            break;
918                                                                    }
919                                                            }
920                                                    }
921                                            }
922                                    }
923                            }
924                            // second pass (ascending the diagonal):
925                            for (int row = (int) source.getRowCount() - 1; row > 0; --row) {
926                                    for (int colData = 0; colData < row - 1; ++colData) {
927                                            if (0 != source.getAsDouble(row, colData)) {
928                                                    for (int colEmpty = colData + 1; colEmpty < row; ++colEmpty) {
929                                                            if (0 == source.getAsDouble(row, colEmpty)) {
930                                                                    if (Ginv.canSwapCols(source, colData, colEmpty, row)) {
931                                                                            Ginv.swapRows(source, colData, colEmpty);
932                                                                            Ginv.swapCols(source, colData, colEmpty);
933                                                                            Ginv.swapRows(response, colData, colEmpty);
934                                                                            break;
935                                                                    }
936                                                            }
937                                                    }
938                                            }
939                                    }
940                            }
941                    }
942                    return response;
943            }
944    
945            /**
946             * Mathematical operator to reduce the bandwidth of a HusoMatrix. The
947             * HusoMatrix must be a square HusoMatrix or no operations are performed.
948             * 
949             * This method reduces a sparse HusoMatrix and returns the numbering of
950             * nodes to achieve this banding. It may be advantageous to run this twice,
951             * though sample cases haven't shown the need. Rows are numbered beginning
952             * with 0. The return HusoMatrix is a vector with what should be used as the
953             * new numbering to achieve minimum banding.
954             * 
955             * Each node in a typical finite-element grid is connected to surrounding
956             * nodes which are connected back to this node. This routine was designed
957             * with this requirement in mind. It may work where a node is connected to
958             * an adjacent node that is not connected back to this node... and this is
959             * quite possible when the grid is on a sphere and the connections are
960             * determined based on initial headings or bearings.
961             * 
962             * @return the vector indicating the numbering required to achieve a minimum
963             *         banding
964             */
965            public static Matrix reduce(Matrix source) {
966                    Matrix response = Matrix.factory.zeros(source.getRowCount(), 1);
967                    for (int row = 0; row < source.getRowCount(); ++row) {
968                            response.setAsDouble(row, row, 0);
969                    }
970                    return source.getRowCount() == source.getColumnCount() ? Ginv.reduce(source, response)
971                                    : response;
972            }
973    
974            /**
975             * Calculate the arbitrariness of the solution. This is a way to find out
976             * how unique the existing inverse is. The equation is here: A * x = b And
977             * the solution is: x = A12 * b + an arbitrariness which could be infinite,
978             * but will follow a general pattern. For instance, if the solution is a
979             * line, it could be anchored in the Y at any arbitrary distance. If the
980             * solution is a plane it could be arbitrarily set to any place in perhaps
981             * two different dimensions.
982             * 
983             * The arbitrariness is calculated by taking the difference between the
984             * complete inverse times the original and subtracting the generalized
985             * inverse times the original matrix. That's the idea, here's the formula:
986             * 
987             * x = A12 * b + (I - (A12 * A)) * z The z is a completely arbitrary vector
988             * (you decide that one). The product (A12 * A) could be the Identity
989             * HusoMatrix, if the solution is unique, in which case the right side drops
990             * out: (I - I) * z
991             * 
992             * Again, it's a lot easier to refer to the http://aktzin.com/ site for the
993             * description and a different way to get this information.
994             * 
995             * @return the matrix (I - (A12 * A))
996             */
997            public static Matrix arbitrariness(Matrix source, Matrix inverse) {
998                    Matrix intermediate = inverse.mtimes(source);
999                    return Matrix.factory.eye(intermediate.getSize()).minus(intermediate);
1000            }
1001    
1002            public double getDouble(long... coordinates) throws MatrixException {
1003                    throw new MatrixException("this method should never be called: LINK not possible");
1004            }
1005    
1006            public DenseDoubleMatrix2D calcLink() throws MatrixException {
1007                    throw new MatrixException("linking not possible, use ORIG or NEW");
1008            }
1009    
1010            public DenseDoubleMatrix2D calcNew() throws MatrixException {
1011                    Matrix source = getSource();
1012                    ArrayDenseDoubleMatrix2D matrix = new ArrayDenseDoubleMatrix2D(source);
1013                    return inverse(matrix.getRowMajorDoubleArray2D());
1014            }
1015    
1016            public DenseDoubleMatrix2D calcOrig() throws MatrixException {
1017                    Matrix source = getSource();
1018                    if (!source.isSquare()) {
1019                            throw new MatrixException("ORIG only possible for square matrices");
1020                    }
1021    
1022                    if (source instanceof HasRowMajorDoubleArray2D) {
1023                            return inverse(((HasRowMajorDoubleArray2D) source).getRowMajorDoubleArray2D());
1024                    } else if (source instanceof DenseDoubleMatrix2D) {
1025                            return inverse((DenseDoubleMatrix2D) source);
1026                    } else {
1027                            return inverse(source);
1028                    }
1029            }
1030    
1031    }