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 }