void purge(int n,int ll,double[] r,double[] q,double[] ra,double[] qa,double[] wrk,double[] eta,double[] oldeta,int step,double[] rnmp,double tol){
  double t, tq, tr, reps1;
  double rnm=rnmp[0];
  int k, iteration, i;
  boolean flag;
  if (step < ll + 2)   return;
  k=svd_idamax(step - (ll + 1),eta,ll,1) + ll;
  if (Math.abs(eta[k]) > reps) {
    reps1=eps1 / reps;
    iteration=0;
    flag=true;
    while (iteration < 2 && flag) {
      if (rnm > tol) {
        tq=0.0;
        tr=0.0;
        for (i=ll; i < step; i++) {
          store(n,RETRQ,i,wrk);
          t=-svd_ddot(n,qa,1,wrk,1);
          tq+=Math.abs(t);
          svd_daxpy(n,t,wrk,1,q,1);
          t=-svd_ddot(n,ra,1,wrk,1);
          tr+=Math.abs(t);
          svd_daxpy(n,t,wrk,1,r,1);
        }
        svd_dcopy(n,q,1,qa,1);
        t=-svd_ddot(n,r,1,qa,1);
        tr+=Math.abs(t);
        svd_daxpy(n,t,q,1,r,1);
        svd_dcopy(n,r,1,ra,1);
        rnm=Math.sqrt(svd_ddot(n,ra,1,r,1));
        if (tq <= reps1 && tr <= reps1 * rnm)         flag=false;
      }
      iteration++;
    }
    for (i=ll; i <= step; i++) {
      eta[i]=eps1;
      oldeta[i]=eps1;
    }
  }
  rnmp[0]=rnm;
  return;
}

void purge ( int n , int ll , double [ ] r , double [ ] q , double [ ] ra , double [ ] qa , double [ ] wrk , double [ ] eta , double [ ] oldeta , int step , double [ ] rnmp , double tol ) { double t , tq , tr , reps1 ; double rnm = rnmp [ 0 ] ; int k , iteration , i ; boolean flag ; if ( step < ll + 2 ) return ; k = svd_idamax ( step - ( ll + 1 ) , eta , ll , 1 ) + ll ; if ( Math . abs ( eta [ k ] ) > reps ) { reps1 = eps1 / reps ; iteration = 0 ; flag = true ; while ( iteration < 2 && flag ) { if ( rnm > tol ) { tq = 0.0 ; tr = 0.0 ; for ( i = ll ; i < step ; i ++ ) { store ( n , RETRQ , i , wrk ) ; t = - svd_ddot ( n , qa , 1 , wrk , 1 ) ; tq += Math . abs ( t ) ; svd_daxpy ( n , t , wrk , 1 , q , 1 ) ; t = - svd_ddot ( n , ra , 1 , wrk , 1 ) ; tr += Math . abs ( t ) ; svd_daxpy ( n , t , wrk , 1 , r , 1 ) ; } svd_dcopy ( n , q , 1 , qa , 1 ) ; t = - svd_ddot ( n , r , 1 , qa , 1 ) ; tr += Math . abs ( t ) ; svd_daxpy ( n , t , q , 1 , r , 1 ) ; svd_dcopy ( n , r , 1 , ra , 1 ) ; rnm = Math . sqrt ( svd_ddot ( n , ra , 1 , r , 1 ) ) ; if ( tq <= reps1 && tr <= reps1 * rnm ) flag = false ; } iteration ++ ; } for ( i = ll ; i <= step ; i ++ ) { eta [ i ] = eps1 ; oldeta [ i ] = eps1 ; } } rnmp [ 0 ] = rnm ; return ; }
description - - - - - - - - - - - function examines the state of orthogonality between the new lanczos vector and the previous ones to decide whether re - orthogonalization should be performed arguments - - - - - - - - - ( input ) n dimension of the eigenproblem for matrix b ll number of intitial lanczos vectors in local orthog . r residual vector to become next lanczos vector q current lanczos vector ra previous lanczos vector qa previous lanczos vector wrk temporary vector to hold the previous lanczos vector eta state of orthogonality between r and prev . lanczos vectors oldeta state of orthogonality between q and prev . lanczos vectors j current lanczos step ( output ) r residual vector orthogonalized against previous lanczos vectors q current lanczos vector orthogonalized against previous ones functions used - - - - - - - - - - - - - - blas svd_daxpy , svd_dcopy , svd_idamax , svd_ddot user store
