Note that the above strategy assumes that the matrix is used for an affine transform, which is not always + * true (it could be the matrix of a map projection derivative for instance). If the matrix is not for an affine + * transform, then the last column has no special meaning and the above strategy is somewhat asymmetric. + * However it will still produce NaN for the full row in matrix multiplications.

+ * + *Conversely, if the matrix has more rows than columns (in a system of linear equations, the system would + * be overdetermined), then we omit the rows containing only zero or NaN values. After the matrix + * inversion, we insert columns having only zero values for the dimensions associated to those rows. + * Semantically, the inverse matrix is a (x₁,y₁,z,t) → (x₂,y₂) transform that just discards the ordinate values + * at the dimensions corresponding to those rows.

*/ @Override public MatrixSIS inverse() throws NoninvertibleMatrixException { + if (numRow < numCol) { + return inverseDimensionReduction(); + } else { + return inverseDimensionIncrease(); + } + } + + /** + * Inverses a matrix for a transform where target points have fewer ordinates than source points. + * If a column contains only zero values, then this means that the ordinate at the corresponding + * column is simply deleted. We can omit that column. We check the last columns before the first + * columns on the assumption that last dimensions are more likely to be independent dimensions + * like time. + */ + private MatrixSIS inverseDimensionReduction() throws NoninvertibleMatrixException { final int numRow = this.numRow; // Protection against accidental changes. final int numCol = this.numCol; - if (numRow < numCol) { - final int length = numRow * numCol; - /* - * Target points have fewer ordinates than source points. If a column contains only zero values, - * then this means that the ordinate at the corresponding column is simply deleted. We can omit - * that column. We check the last columns before the first columns on the assumption that last - * dimensions are more likely to be independant dimensions like time. - */ - int oi = numCol - numRow; - final int[] omitted = new int[oi]; -skipColumn: for (int i=numCol; --i>=0;) { - for (int j=length + i; (j -= numCol) >= 0;) { - if (elements[j] != 0) { - continue skipColumn; - } + final int length = numRow * numCol; + int i = numCol; + int oi = numCol - numRow; // Initialized to the maximal amount of columns that we may omit. + final int[] omitted = new int[oi]; +next: do { + if (--i < 0) { + throw nonInvertible(); // Not enough columns that we can omit. + } + for (int j=length + i; (j -= numCol) >= 0;) { + if (elements[j] != 0) { + continue next; } - omitted[--oi] = i; // Found a column which contains only 0 elements. - if (oi == 0) { - /* - * Found enough columns containing only zero elements. Create a square matrix omitting those - * columns, and invert that matrix. Note that we also need to either copy the error terms, - * or to infer them. - */ - GeneralMatrix squareMatrix = new GeneralMatrix(numRow, numRow, false, 2); - int j=0; - for (i=0; i