/[escript]/branches/doubleplusgood/paso/src/TFQMR.cpp
ViewVC logotype

Diff of /branches/doubleplusgood/paso/src/TFQMR.cpp

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

trunk/paso/src/TFQMR.c revision 2826 by artak, Fri Dec 18 01:33:35 2009 UTC branches/doubleplusgood/paso/src/TFQMR.cpp revision 4275 by jfenwick, Tue Mar 5 09:32:52 2013 UTC
# Line 1  Line 1 
1    
2  /*******************************************************  /*****************************************************************************
3  *  *
4  * Copyright (c) 2003-2009 by University of Queensland  * Copyright (c) 2003-2013 by University of Queensland
5  * Earth Systems Science Computational Center (ESSCC)  * http://www.uq.edu.au
 * http://www.uq.edu.au/esscc  
6  *  *
7  * Primary Business: Queensland, Australia  * Primary Business: Queensland, Australia
8  * Licensed under the Open Software License version 3.0  * Licensed under the Open Software License version 3.0
9  * http://www.opensource.org/licenses/osl-3.0.php  * http://www.opensource.org/licenses/osl-3.0.php
10  *  *
11  *******************************************************/  * Development until 2012 by Earth Systems Science Computational Center (ESSCC)
12    * Development since 2012 by School of Earth Sciences
13    *
14    *****************************************************************************/
15    
16    
17  /* TFQMR iterations */  /* TFQMR iterations */
# Line 23  Line 25 
25  #include <omp.h>  #include <omp.h>
26  #endif  #endif
27    
28  #ifdef PASO_MPI  #ifdef ESYS_MPI
29  #include <mpi.h>  #include <mpi.h>
30  #endif  #endif
31    
# Line 53  Line 55 
55  *  INFO    (output) INT  *  INFO    (output) INT
56  *  *
57  *          = SOLVER_NO_ERROR: Successful exit. Iterated approximate solution returned.  *          = SOLVER_NO_ERROR: Successful exit. Iterated approximate solution returned.
58  *          = SOLVEr_MAXITER_REACHED  *          = SOLVER_MAXITER_REACHED
59  *          = SOLVER_INPUT_ERROR Illegal parameter:  *          = SOLVER_INPUT_ERROR Illegal parameter:
60  *          = SOLVEr_BREAKDOWN: If parameters rHO or OMEGA become smaller  *          = SOLVER_BREAKDOWN: If parameters RHO or OMEGA become smaller
61  *          = SOLVER_MEMORY_ERROR : If parameters rHO or OMEGA become smaller  *          = SOLVER_MEMORY_ERROR : If parameters RHO or OMEGA become smaller
62  *  *
63  *  ==============================================================  *  ==============================================================
64  */  */
# Line 84  err_t Paso_Solver_TFQMR( Line 86  err_t Paso_Solver_TFQMR(
86    bool_t breakFlag=FALSE, maxIterFlag=FALSE, convergeFlag=FALSE;    bool_t breakFlag=FALSE, maxIterFlag=FALSE, convergeFlag=FALSE;
87    err_t status = SOLVER_NO_ERROR;    err_t status = SOLVER_NO_ERROR;
88    dim_t n = Paso_SystemMatrix_getTotalNumRows(A);    dim_t n = Paso_SystemMatrix_getTotalNumRows(A);
89    double  *u1=NULL, *u2=NULL, *y1=NULL, *y2=NULL, *d=NULL, *w=NULL, *v=NULL, *v_old=NULL, *temp_vector=NULL,*res=NULL;    double  *u1=NULL, *u2=NULL, *y1=NULL, *y2=NULL, *d=NULL, *w=NULL, *v=NULL, *temp_vector=NULL,*res=NULL;
90    
91    double eta,theta,tau,rho,beta,alpha,sigma,rhon,c;    double eta,theta,tau,rho,beta,alpha,sigma,rhon,c;
92    
# Line 97  err_t Paso_Solver_TFQMR( Line 99  err_t Paso_Solver_TFQMR(
99  /*   ---------------------                                         */  /*   ---------------------                                         */
100  /*                                                                 */  /*                                                                 */
101  /*                                                                 */  /*                                                                 */
102    u1=TMPMEMALLOC(n,double);    u1=new double[n];
103    u2=TMPMEMALLOC(n,double);    u2=new double[n];
104    y1=TMPMEMALLOC(n,double);    y1=new double[n];
105    y2=TMPMEMALLOC(n,double);    y2=new double[n];
106    d=TMPMEMALLOC(n,double);    d=new double[n];
107    w=TMPMEMALLOC(n,double);    w=new double[n];
108    v=TMPMEMALLOC(n,double);    v=new double[n];
109    v_old=TMPMEMALLOC(n,double);    temp_vector=new double[n];
110    temp_vector=TMPMEMALLOC(n,double);    res=new double[n];
   res=TMPMEMALLOC(n,double);  
111    
112   if (u1 ==NULL || u2== NULL || y1 == NULL || y2== NULL || d==NULL || w==NULL || v==NULL || v_old==NULL) {   if (u1 ==NULL || u2== NULL || y1 == NULL || y2== NULL || d==NULL || w==NULL || v==NULL ) {
113       status=SOLVER_MEMORY_ERROR;       status=SOLVER_MEMORY_ERROR;
114    }    }
115    
# Line 122  err_t Paso_Solver_TFQMR( Line 123  err_t Paso_Solver_TFQMR(
123    Paso_zeroes(n,x);    Paso_zeroes(n,x);
124        
125    Performance_startMonitor(pp,PERFORMANCE_PRECONDITIONER);    Performance_startMonitor(pp,PERFORMANCE_PRECONDITIONER);
126    Paso_Solver_solvePreconditioner(A,res,r);    Paso_SystemMatrix_solvePreconditioner(A,res,r);
127    Performance_stopMonitor(pp,PERFORMANCE_PRECONDITIONER);    Performance_stopMonitor(pp,PERFORMANCE_PRECONDITIONER);
128        
129    Performance_startMonitor(pp,PERFORMANCE_SOLVER);    Performance_startMonitor(pp,PERFORMANCE_SOLVER);
# Line 143  err_t Paso_Solver_TFQMR( Line 144  err_t Paso_Solver_TFQMR(
144        
145    Performance_stopMonitor(pp,PERFORMANCE_SOLVER);    Performance_stopMonitor(pp,PERFORMANCE_SOLVER);
146    Performance_startMonitor(pp,PERFORMANCE_PRECONDITIONER);    Performance_startMonitor(pp,PERFORMANCE_PRECONDITIONER);
147    Paso_Solver_solvePreconditioner(A,v,temp_vector);    Paso_SystemMatrix_solvePreconditioner(A,v,temp_vector);
148    Performance_stopMonitor(pp,PERFORMANCE_PRECONDITIONER);    Performance_stopMonitor(pp,PERFORMANCE_PRECONDITIONER);
149    Performance_startMonitor(pp,PERFORMANCE_SOLVER);    Performance_startMonitor(pp,PERFORMANCE_SOLVER);
150      /* v = P^{-1} * A y1 */
151        
152    Paso_Copy(n,u1,v);    Paso_Copy(n,u1,v);
153            
# Line 156  err_t Paso_Solver_TFQMR( Line 158  err_t Paso_Solver_TFQMR(
158        
159    rho = tau * tau;    rho = tau * tau;
160                
161    norm_of_residual=tau*sqrt ( m + 1 );    norm_of_residual=tau;
162        
163    while (!(convergeFlag || maxIterFlag || breakFlag || (status !=SOLVER_NO_ERROR) ))    while (!(convergeFlag || maxIterFlag || breakFlag || (status !=SOLVER_NO_ERROR) ))
164    {    {
# Line 172  err_t Paso_Solver_TFQMR( Line 174  err_t Paso_Solver_TFQMR(
174         {         {
175           /*  Compute y2 and u2 only if you have to */           /*  Compute y2 and u2 only if you have to */
176           if ( j == 1 ){           if ( j == 1 ){
177            Paso_LinearCombination(n,y2,1.,y1,-alpha,v);            Paso_LinearCombination(n,y2,PASO_ONE,y1,-alpha,v); /* y2 = y1 - alpha*v */
178                        
179            Performance_stopMonitor(pp,PERFORMANCE_SOLVER);            Performance_stopMonitor(pp,PERFORMANCE_SOLVER);
180            Performance_startMonitor(pp,PERFORMANCE_MVM);            Performance_startMonitor(pp,PERFORMANCE_MVM);
# Line 182  err_t Paso_Solver_TFQMR( Line 184  err_t Paso_Solver_TFQMR(
184                        
185            Performance_stopMonitor(pp,PERFORMANCE_SOLVER);            Performance_stopMonitor(pp,PERFORMANCE_SOLVER);
186            Performance_startMonitor(pp,PERFORMANCE_PRECONDITIONER);            Performance_startMonitor(pp,PERFORMANCE_PRECONDITIONER);
187            Paso_Solver_solvePreconditioner(A,u2,temp_vector);        Paso_SystemMatrix_solvePreconditioner(A,u2,temp_vector);  
188            Performance_stopMonitor(pp,PERFORMANCE_PRECONDITIONER);            Performance_stopMonitor(pp,PERFORMANCE_PRECONDITIONER);
189            Performance_startMonitor(pp,PERFORMANCE_SOLVER);            Performance_startMonitor(pp,PERFORMANCE_SOLVER);
190              /* u2 = P^{-1} * A y2 */
191           }           }
192           m = 2 * (num_iter+1) - 2 + (j+1);           m = 2 * (num_iter+1) - 2 + (j+1);
193    
194            if (j==0) {            if (j==0) {
195              Paso_Update(n,1.,w,-alpha,u1);              Paso_Update(n,1.,w,-alpha,u1); /* w = w - alpha * u1 */
196              Paso_Update(n,( theta * theta * eta / alpha ),d,1.,y1);              Paso_Update(n,( theta * theta * eta / alpha ),d,1.,y1); /* d = ( theta * theta * eta / alpha )*d + y1 */
197            }            }
198            if (j==1) {            if (j==1) {
199              Paso_Update(n,1.,w,-alpha,u2);              Paso_Update(n,1.,w,-alpha,u2);  /* w = w - -alpha * u2 */
200              Paso_Update(n,( theta * theta * eta / alpha ),d,1.,y2);              Paso_Update(n,( theta * theta * eta / alpha ),d,1.,y2); /* d = ( theta * theta * eta / alpha )*d + y2 */
201            }            }
202                    
203           theta =Paso_l2(n,w,A->mpi_info)/tau;           theta =Paso_l2(n,w,A->mpi_info)/tau;
204           c = 1.0 / sqrt ( 1.0 + theta * theta );           /*printf("tau = %e, %e %e\n",tau, Paso_l2(n,w,A->mpi_info)/tau, theta);*/
205             c = PASO_ONE / sqrt ( PASO_ONE + theta * theta );
206           tau = tau * theta * c;           tau = tau * theta * c;
207           eta = c * c * alpha;           eta = c * c * alpha;
208           Paso_Update(n,1.,x,eta,d);                 Paso_Update(n,1.,x,eta,d);      
# Line 209  err_t Paso_Solver_TFQMR( Line 214  err_t Paso_Solver_TFQMR(
214       beta = rhon / rho;       beta = rhon / rho;
215       rho = rhon;       rho = rhon;
216        
217       Paso_LinearCombination(n,y1,1.,w,beta,y2);       Paso_LinearCombination(n,y1, PASO_ONE,w,beta,y2); /* y1 = w + beta * y2 */
218    
219       Performance_stopMonitor(pp,PERFORMANCE_SOLVER);       Performance_stopMonitor(pp,PERFORMANCE_SOLVER);
       
220       Performance_startMonitor(pp,PERFORMANCE_MVM);       Performance_startMonitor(pp,PERFORMANCE_MVM);
221       Paso_SystemMatrix_MatrixVector_CSR_OFFSET0(PASO_ONE, A, y1,PASO_ZERO,temp_vector);       Paso_SystemMatrix_MatrixVector_CSR_OFFSET0(PASO_ONE, A, y1,PASO_ZERO,temp_vector);
222       Performance_stopMonitor(pp,PERFORMANCE_MVM);       Performance_stopMonitor(pp,PERFORMANCE_MVM);
223            
224       Performance_startMonitor(pp,PERFORMANCE_PRECONDITIONER);       Performance_startMonitor(pp,PERFORMANCE_PRECONDITIONER);
225       Paso_Solver_solvePreconditioner(A,u1,temp_vector);       Paso_SystemMatrix_solvePreconditioner(A,u1,temp_vector);
226       Performance_stopMonitor(pp,PERFORMANCE_PRECONDITIONER);       Performance_stopMonitor(pp,PERFORMANCE_PRECONDITIONER);
       
227       Performance_startMonitor(pp,PERFORMANCE_SOLVER);       Performance_startMonitor(pp,PERFORMANCE_SOLVER);
228         /*  u1 = P^{-1} * A y1 */
229            
230       Paso_Copy(n,v_old,v);       Paso_LinearCombination(n,temp_vector,PASO_ONE,u2,beta,v); /* t = u2 + beta * v */
231                   Paso_LinearCombination(n,v,PASO_ONE,u1,beta,temp_vector); /* v = u1 + beta * t */
      Paso_Update(n,beta,v_old,1,u2);  
      Paso_LinearCombination(n,v,1.,u1,beta,v_old);  
232       }       }
233       maxIterFlag = (num_iter > maxit);       maxIterFlag = (num_iter > maxit);
234       norm_of_residual=tau*sqrt ( m + 1 );       norm_of_residual=tau*sqrt ( (double) (m + 1 ) );
235       convergeFlag=(norm_of_residual<(*tolerance));       convergeFlag=(norm_of_residual<(*tolerance));
236            
237            
# Line 243  err_t Paso_Solver_TFQMR( Line 245  err_t Paso_Solver_TFQMR(
245      /* end of iteration */      /* end of iteration */
246            
247      Performance_stopMonitor(pp,PERFORMANCE_SOLVER);      Performance_stopMonitor(pp,PERFORMANCE_SOLVER);
248      TMPMEMFREE(u1);      delete[] (u1);
249      TMPMEMFREE(u2);      delete[] (u2);
250      TMPMEMFREE(y1);      delete[] (y1);
251      TMPMEMFREE(y2);      delete[] (y2);
252      TMPMEMFREE(d);      delete[] (d);
253      TMPMEMFREE(w);      delete[] (w);
254      TMPMEMFREE(v);      delete[] (v);
255      TMPMEMFREE(v_old);      delete[] (temp_vector);
256      TMPMEMFREE(temp_vector);      delete[] (res);
     TMPMEMFREE(res);  
257      *iter=num_iter;      *iter=num_iter;
258      *tolerance=norm_of_residual;      *tolerance=norm_of_residual;
259            

Legend:
Removed from v.2826  
changed lines
  Added in v.4275

  ViewVC Help
Powered by ViewVC 1.1.26