/[escript]/trunk/paso/src/MINRES.c
ViewVC logotype

Diff of /trunk/paso/src/MINRES.c

Parent Directory Parent Directory | Revision Log Revision Log | View Patch Patch

revision 1787 by artak, Mon Sep 15 01:36:34 2008 UTC revision 2565 by artak, Tue Jul 28 05:50:15 2009 UTC
# Line 1  Line 1 
1    
2  /*******************************************************  /*******************************************************
3   *  *
4   *           Copyright 2003-2007 by ACceSS MNRF  * Copyright (c) 2003-2009 by University of Queensland
5   *       Copyright 2007 by University of Queensland  * Earth Systems Science Computational Center (ESSCC)
6   *  * http://www.uq.edu.au/esscc
7   *                http://esscc.uq.edu.au  *
8   *        Primary Business: Queensland, Australia  * Primary Business: Queensland, Australia
9   *  Licensed under the Open Software License version 3.0  * Licensed under the Open Software License version 3.0
10   *     http://www.opensource.org/licenses/osl-3.0.php  * http://www.opensource.org/licenses/osl-3.0.php
11   *  *
12   *******************************************************/  *******************************************************/
13    
14  /* TFQMR iterations */  
15    /* MINRES iterations */
16    
17  #include "SystemMatrix.h"  #include "SystemMatrix.h"
18  #include "Paso.h"  #include "Paso.h"
# Line 71  err_t Paso_Solver_MINRES( Line 72  err_t Paso_Solver_MINRES(
72      double * r,      double * r,
73      double * x,      double * x,
74      dim_t *iter,      dim_t *iter,
     double *tol,  
75      double * tolerance,      double * tolerance,
76      Paso_Performance* pp) {      Paso_Performance* pp) {
77    
# Line 83  err_t Paso_Solver_MINRES( Line 83  err_t Paso_Solver_MINRES(
83    dim_t n = Paso_SystemMatrix_getTotalNumRows(A);    dim_t n = Paso_SystemMatrix_getTotalNumRows(A);
84    double  *w=NULL, *w1=NULL, *w2=NULL, *r1=NULL, *r2=NULL, *y=NULL, *v=NULL;    double  *w=NULL, *w1=NULL, *w2=NULL, *r1=NULL, *r2=NULL, *y=NULL, *v=NULL;
85    
86    double Anorm,ynorm,oldb,dbar,epsln,phibar,rhs1,rhs2,rnorm,tnorm2,ynorm2,cs,sn,eps,s,alfa,denom,z,beta1,beta;    double Anorm,Arnorm,ynorm,oldb,dbar,epsln,phibar,rhs1,rhs2,rnorm,tnorm2,ynorm2,cs,sn,eps,s,alfa,denom,z,beta1,beta;
87    double gmax,gmin,oldeps,delta,gbar,gamma,phi;    double gmax,gmin,oldeps,delta,gbar,gamma,phi,root,epsx;
88    
89    double norm_of_residual;    double norm_of_residual=0;
90        
91  /*                                                                 */  /*                                                                 */
92  /*-----------------------------------------------------------------*/  /*-----------------------------------------------------------------*/
# Line 151  err_t Paso_Solver_MINRES( Line 151  err_t Paso_Solver_MINRES(
151    ynorm2 = 0;    ynorm2 = 0;
152    cs     = -1;    cs     = -1;
153    sn     = 0;    sn     = 0;
154    eps    = 0.0001;    eps    = 0.000001;
155    
156    while (!(convergeFlag || maxIterFlag || breakFlag || (status !=SOLVER_NO_ERROR) ))    while (!(convergeFlag || maxIterFlag || breakFlag || (status !=SOLVER_NO_ERROR) ))
157    {    {
# Line 161  err_t Paso_Solver_MINRES( Line 161  err_t Paso_Solver_MINRES(
161            
162       Performance_stopMonitor(pp,PERFORMANCE_SOLVER);       Performance_stopMonitor(pp,PERFORMANCE_SOLVER);
163       Performance_startMonitor(pp,PERFORMANCE_MVM);       Performance_startMonitor(pp,PERFORMANCE_MVM);
164       Paso_SystemMatrix_MatrixVector_CSR_OFFSET0(ONE, A, v,ZERO,y);       Paso_SystemMatrix_MatrixVector_CSR_OFFSET0(PASO_ONE, A, v,PASO_ZERO,y);
165       Performance_stopMonitor(pp,PERFORMANCE_MVM);       Performance_stopMonitor(pp,PERFORMANCE_MVM);
166       Performance_startMonitor(pp,PERFORMANCE_SOLVER);       Performance_startMonitor(pp,PERFORMANCE_SOLVER);
167            
# Line 170  err_t Paso_Solver_MINRES( Line 170  err_t Paso_Solver_MINRES(
170       }       }
171    
172       alfa = Paso_InnerProduct(n,v,y,A->mpi_info);       alfa = Paso_InnerProduct(n,v,y,A->mpi_info);
173       Paso_Update(n, 1., y, -(alfa/beta), r2);       Paso_Update(n, 1., y, (-alfa/beta), r2);
174       Paso_Copy(n,r1,r2);       Paso_Copy(n,r1,r2);
175       Paso_Copy(n,r2,y);       Paso_Copy(n,r2,y);
176    
# Line 204  err_t Paso_Solver_MINRES( Line 204  err_t Paso_Solver_MINRES(
204       epsln  =               sn * beta ;       epsln  =               sn * beta ;
205       dbar   =            -  cs * beta;       dbar   =            -  cs * beta;
206            
207         root   = sqrt(gbar*gbar+dbar*dbar) ;
208         Arnorm = phibar*root;
209        
210       gamma  = sqrt(gbar*gbar+beta*beta) ;       gamma  = sqrt(gbar*gbar+beta*beta) ;
211       gamma  = MAX(gamma,eps) ;       gamma  = MAX(gamma,eps) ;
212       cs     = gbar / gamma ;                   cs     = gbar / gamma ;            
# Line 234  err_t Paso_Solver_MINRES( Line 237  err_t Paso_Solver_MINRES(
237       ynorm  = sqrt( ynorm2 ) ;       ynorm  = sqrt( ynorm2 ) ;
238    
239       rnorm  = phibar;       rnorm  = phibar;
240         epsx   = Anorm*ynorm*eps;
241    
242       maxIterFlag = (num_iter > maxit);       maxIterFlag = (num_iter > maxit);
243       norm_of_residual=rnorm;       norm_of_residual=rnorm;
244       convergeFlag=(norm_of_residual<Anorm*ynorm*(*tolerance));       convergeFlag=((norm_of_residual/(Anorm*ynorm))<(*tolerance) || 1+(norm_of_residual/(Anorm*ynorm)) <=1);
       
245            
246       if (maxIterFlag) {       if (maxIterFlag) {
247           status = SOLVER_MAXITER_REACHED;           status = SOLVER_MAXITER_REACHED;
# Line 246  err_t Paso_Solver_MINRES( Line 249  err_t Paso_Solver_MINRES(
249           status = SOLVER_BREAKDOWN;           status = SOLVER_BREAKDOWN;
250       }       }
251      ++(num_iter);      ++(num_iter);
     /*printf("residual norm %.10f < %.10f %.10f %.10f \n",rnorm,Anorm*ynorm*(*tolerance), Anorm*ynorm, (*tolerance));*/  
252    }    }
253      /* end of iteration */      /* end of iteration */
254            
# Line 260  err_t Paso_Solver_MINRES( Line 262  err_t Paso_Solver_MINRES(
262      TMPMEMFREE(v);      TMPMEMFREE(v);
263        
264      *iter=num_iter;      *iter=num_iter;
265      *tol=norm_of_residual;      *tolerance=norm_of_residual;
266            
267    /*     End of MINRES */    /*     End of MINRES */
268    return status;    return status;

Legend:
Removed from v.1787  
changed lines
  Added in v.2565

  ViewVC Help
Powered by ViewVC 1.1.26