@ Override public DoubleMatrix1D [ ] getMatrixScalingFactors ( DoubleMatrix2D A ) { DoubleFactory1D F1 = DoubleFactory1D . dense ; Algebra ALG = Algebra . DEFAULT ; int r = A . rows ( ) ; int c = A . columns ( ) ; DoubleMatrix1D D1 = F1 . make ( r , 1 ) ; DoubleMatrix1D D2 = F1 . make ( c , 1 ) ; DoubleMatrix2D AK = A . copy ( ) ; DoubleMatrix1D DR = F1 . make ( r , 1 ) ; DoubleMatrix1D DC = F1 . make ( c , 1 ) ; DoubleMatrix1D DRInv = F1 . make ( r ) ; DoubleMatrix1D DCInv = F1 . make ( c ) ; log . debug ( "eps : " + eps ) ; int maxIteration = 50 ; for ( int k = 0 ; k <= maxIteration ; k ++ ) { double normR = - Double . MAX_VALUE ; double normC = - Double . MAX_VALUE ; for ( int i = 0 ; i < r ; i ++ ) { double dri = ALG . normInfinity ( AK . viewRow ( i ) ) ; DR . setQuick ( i , Math . sqrt ( dri ) ) ; DRInv . setQuick ( i , 1. / Math . sqrt ( dri ) ) ; normR = Math . max ( normR , Math . abs ( 1 - dri ) ) ; } for ( int j = 0 ; j < c ; j ++ ) { double dci = ALG . normInfinity ( AK . viewColumn ( j ) ) ; DC . setQuick ( j , Math . sqrt ( dci ) ) ; DCInv . setQuick ( j , 1. / Math . sqrt ( dci ) ) ; normC = Math . max ( normC , Math . abs ( 1 - dci ) ) ; } log . debug ( "normR: " + normR ) ; log . debug ( "normC: " + normC ) ; if ( normR < eps && normC < eps ) { break ; } //D1 = ALG.mult(D1, DRInv); for ( int i = 0 ; i < r ; i ++ ) { double prevD1I = D1 . getQuick ( i ) ; double newD1I = prevD1I * DRInv . getQuick ( i ) ; D1 . setQuick ( i , newD1I ) ; } //D2 = ALG.mult(D2, DCInv); for ( int j = 0 ; j < c ; j ++ ) { double prevD2J = D2 . getQuick ( j ) ; double newD2J = prevD2J * DCInv . getQuick ( j ) ; D2 . setQuick ( j , newD2J ) ; } //log.debug("D1: " + ArrayUtils.toString(D1.toArray())); //log.debug("D2: " + ArrayUtils.toString(D2.toArray())); if ( k == maxIteration ) { log . warn ( "max iteration reached" ) ; } //AK = ALG.mult(DRInv, ALG.mult(AK, DCInv)); AK = ColtUtils . diagonalMatrixMult ( DRInv , AK , DCInv ) ; } return new DoubleMatrix1D [ ] { D1 , D2 } ; }
Scaling factors for not singular matrices .