public QRP factor ( Matrix A ) { if ( Q . numRows ( ) != A . numRows ( ) ) throw new IllegalArgumentException ( "Q.numRows() != A.numRows()" ) ; else if ( R . numColumns ( ) != A . numColumns ( ) ) throw new IllegalArgumentException ( "R.numColumns() != A.numColumns()" ) ; // copy A values in Afact Afact . zero ( ) ; for ( MatrixEntry e : A ) { Afact . set ( e . row ( ) , e . column ( ) , e . get ( ) ) ; } intW info = new intW ( 0 ) ; LAPACK lapack = LAPACK . getInstance ( ) ; // the dgeqp3 work array is a different size than the dorgqr // one (I'm not entirely convinced that we should be // persisting the work arrays as fields). double [ ] factorWorkOptimalSize = { 0.0 } ; double [ ] factorWork ; lapack . dgeqp3 ( m , n , Afact . getData ( ) , Matrices . ld ( m ) , jpvt , tau , factorWorkOptimalSize , - 1 , info ) ; factorWork = new double [ ( int ) factorWorkOptimalSize [ 0 ] ] ; /* * Calculate factorisation */ lapack . dgeqp3 ( m , n , Afact . getData ( ) , Matrices . ld ( m ) , jpvt , tau , factorWork , factorWork . length , info ) ; if ( info . val < 0 ) { throw new IllegalArgumentException ( "DGEQP3 was " + info . val ) ; } /* * Get R from Afact */ R . zero ( ) ; for ( MatrixEntry e : Afact ) { if ( e . row ( ) <= e . column ( ) && e . column ( ) < R . numColumns ( ) ) { R . set ( e . row ( ) , e . column ( ) , e . get ( ) ) ; } } /* * Calculate the rank based on a precision EPS */ final double EPS = 1e-12 ; for ( rank = 0 ; rank < k ; rank ++ ) { if ( Math . abs ( R . get ( rank , rank ) ) < EPS ) break ; } /* * Explicit the orthogonal matrix */ lapack . dorgqr ( m , m , k , Afact . getData ( ) , Matrices . ld ( m ) , tau , work , work . length , info ) ; for ( MatrixEntry e : Afact ) { if ( e . column ( ) < Q . numColumns ( ) ) Q . set ( e . row ( ) , e . column ( ) , e . get ( ) ) ; } if ( info . val < 0 ) throw new IllegalArgumentException ( ) ; // Adjust the permutation to zero offset for ( int i = 0 ; i < jpvt . length ; i ++ ) { -- jpvt [ i ] ; } return this ; }
Executes a QR factorization for the given matrix .