private static double purge(int ll, double[][] Q, double[] r, double[] q, double[] ra, double[] qa, double[] eta, double[] oldeta, int step, double rnm, double tol, double eps, double reps) {
        if (step < ll + 2) {
            return rnm;
        }

        double t, tq, tr;
        int k = idamax(step - (ll + 1), eta, ll, 1) + ll;
        if (Math.abs(eta[k]) > reps) {
            double reps1 = eps / reps;
            int iteration = 0;
            boolean flag = true;
            while (iteration < 2 && flag) {
                if (rnm > tol) {
                    // bring in a lanczos vector t and orthogonalize both
                    // r and q against it
                    tq = 0.0;
                    tr = 0.0;
                    for (int i = ll; i < step; i++) {
                        t = -Math.dot(qa, Q[i]);
                        tq += Math.abs(t);
                        Math.axpy(t, Q[i], q);
                        t = -Math.dot(ra, Q[i]);
                        tr += Math.abs(t);
                        Math.axpy(t, Q[i], r);
                    }
                    Math.copy(q, qa);
                    t = -Math.dot(r, qa);
                    tr += Math.abs(t);
                    Math.axpy(t, q, r);
                    Math.copy(r, ra);
                    rnm = Math.sqrt(Math.dot(ra, r));
                    if (tq <= reps1 && tr <= reps1 * rnm) {
                        flag = false;
                    }
                }
                iteration++;
            }

            for (int i = ll; i <= step; i++) {
                eta[i] = eps;
                oldeta[i] = eps;
            }
        }

        return rnm;
    }
private static double purge ( int ll , double [ ] [ ] Q , double [ ] r , double [ ] q , double [ ] ra , double [ ] qa , double [ ] eta , double [ ] oldeta , int step , double rnm , double tol , double eps , double reps ) { if ( step < ll + 2 ) { return rnm ; } double t , tq , tr ; int k = idamax ( step - ( ll + 1 ) , eta , ll , 1 ) + ll ; if ( Math . abs ( eta [ k ] ) > reps ) { double reps1 = eps / reps ; int iteration = 0 ; boolean flag = true ; while ( iteration < 2 && flag ) { if ( rnm > tol ) { tq = 0.0 ; tr = 0.0 ; for ( int i = ll ; i < step ; i ++ ) { t = - Math . dot ( qa , Q [ i ] ) ; tq += Math . abs ( t ) ; Math . axpy ( t , Q [ i ] , q ) ; t = - Math . dot ( ra , Q [ i ] ) ; tr += Math . abs ( t ) ; Math . axpy ( t , Q [ i ] , r ) ; } Math . copy ( q , qa ) ; t = - Math . dot ( r , qa ) ; tr += Math . abs ( t ) ; Math . axpy ( t , q , r ) ; Math . copy ( r , ra ) ; rnm = Math . sqrt ( Math . dot ( ra , r ) ) ; if ( tq <= reps1 && tr <= reps1 * rnm ) { flag = false ; } } iteration ++ ; } for ( int i = ll ; i <= step ; i ++ ) { eta [ i ] = eps ; oldeta [ i ] = eps ; } } return rnm ; }
Examine the state of orthogonality between the new Lanczos
vector and the previous ones to decide whether re-orthogonalization
should be performed.
@param ll       number of intitial Lanczos vectors in local orthog.
@param r        on input, residual vector to become next Lanczos vector.
On output, residual vector orthogonalized against previous Lanczos.
@param q        on input, current Lanczos vector. On Output, current
Lanczos vector orthogonalized against previous ones.
@param ra       previous Lanczos vector
@param qa       previous Lanczos vector
@param eta      state of orthogonality between r and prev. Lanczos vectors
@param oldeta   state of orthogonality between q and prev. Lanczos vectors