public static void invert ( final int blockLength , final boolean upper , final DSubmatrixD1 T , final DSubmatrixD1 T_inv , final double temp [ ] ) { if ( upper ) throw new IllegalArgumentException ( "Upper triangular matrices not supported yet" ) ; if ( temp . length < blockLength * blockLength ) throw new IllegalArgumentException ( "Temp must be at least blockLength*blockLength long." ) ; if ( T . row0 != T_inv . row0 || T . row1 != T_inv . row1 || T . col0 != T_inv . col0 || T . col1 != T_inv . col1 ) throw new IllegalArgumentException ( "T and T_inv must be at the same elements in the matrix" ) ; final int M = T . row1 - T . row0 ; final double dataT [ ] = T . original . data ; final double dataX [ ] = T_inv . original . data ; final int offsetT = T . row0 * T . original . numCols + M * T . col0 ; for ( int i = 0 ; i < M ; i += blockLength ) { int heightT = Math . min ( T . row1 - ( i + T . row0 ) , blockLength ) ; int indexII = offsetT + T . original . numCols * ( i + T . row0 ) + heightT * ( i + T . col0 ) ; for ( int j = 0 ; j < i ; j += blockLength ) { int widthX = Math . min ( T . col1 - ( j + T . col0 ) , blockLength ) ; for ( int w = 0 ; w < temp . length ; w ++ ) { temp [ w ] = 0 ; } for ( int k = j ; k < i ; k += blockLength ) { int widthT = Math . min ( T . col1 - ( k + T . col0 ) , blockLength ) ; int indexL = offsetT + T . original . numCols * ( i + T . row0 ) + heightT * ( k + T . col0 ) ; int indexX = offsetT + T . original . numCols * ( k + T . row0 ) + widthT * ( j + T . col0 ) ; blockMultMinus ( dataT , dataX , temp , indexL , indexX , 0 , heightT , widthT , widthX ) ; } int indexX = offsetT + T . original . numCols * ( i + T . row0 ) + heightT * ( j + T . col0 ) ; InnerTriangularSolver_DDRB . solveL ( dataT , temp , heightT , widthX , heightT , indexII , 0 ) ; System . arraycopy ( temp , 0 , dataX , indexX , widthX * heightT ) ; } InnerTriangularSolver_DDRB . invertLower ( dataT , dataX , heightT , indexII , indexII ) ; } }
Inverts an upper or lower triangular block submatrix .