# Contents of /trunk/paso/src/GMRES.c

Revision 1556 - (show annotations)
Mon May 12 00:54:58 2008 UTC (11 years, 9 months ago) by gross
File MIME type: text/plain
File size: 22221 byte(s)
```Modification to allow mixed mode execution.
In order to keep the code portable accross platform all MPI calls within
parallel regions have been moved.

```
 1 2 /* \$Id\$ */ 3 4 /******************************************************* 5 * 6 * Copyright 2003-2007 by ACceSS MNRF 7 * Copyright 2007 by University of Queensland 8 * 9 * http://esscc.uq.edu.au 10 * Primary Business: Queensland, Australia 11 * Licensed under the Open Software License version 3.0 12 * http://www.opensource.org/licenses/osl-3.0.php 13 * 14 *******************************************************/ 15 16 /* 17 * Purpose 18 * ======= 19 * 20 * GMRES solves the linear system A*x=b using the 21 * truncated and restered GMRES method with preconditioning. 22 * 23 * Convergence test: norm( b - A*x )< TOL. 24 * 25 * 26 * 27 * Arguments 28 * ========= 29 * 30 * r (input/output) double array, dimension n. 31 * On entry, residual of inital guess X 32 * 33 * x (input/output) double array, dimension n. 34 * On input, the initial guess. 35 * 36 * iter (input/output) int 37 * On input, the maximum num_iterations to be performed. 38 * On output, actual number of num_iterations performed. 39 * 40 * tolerance (input/output) DOUBLE PRECISION 41 * On input, the allowable convergence measure for 42 * norm( b - A*x ) 43 * On output, the final value of this measure. 44 * 45 * Length_of_recursion (input) gives the number of residual to be kept in orthogonalization process 46 * 47 * restart (input) If restart>0, iteration is resterted a after restart steps. 48 * 49 * INFO (output) int 50 * 51 * =SOLVER_NO_ERROR: Successful exit. num_iterated approximate solution returned. 52 * =SOLVER_MAXNUM_ITER_REACHED 53 * =SOLVER_INPUT_ERROR Illegal parameter: 54 * =SOLVER_BREAKDOWN: bad luck! 55 * =SOLVER_MEMORY_ERROR : no memory available 56 * 57 * ============================================================== 58 */ 59 60 #include "Common.h" 61 #include "SystemMatrix.h" 62 #include "Solver.h" 63 #ifdef _OPENMP 64 #include 65 #endif 66 67 err_t Paso_Solver_GMRES( 68 Paso_SystemMatrix * A, 69 double * r, 70 double * x, 71 dim_t *iter, 72 double * tolerance,dim_t Length_of_recursion,dim_t restart, 73 Paso_Performance* pp) { 74 75 /* Local variables */ 76 77 double *AP,**X_PRES,**R_PRES,**P_PRES; 78 double *P_PRES_dot_AP,*R_PRES_dot_P_PRES,*BREAKF,*ALPHA; 79 double P_PRES_dot_AP0,P_PRES_dot_AP1,P_PRES_dot_AP2,P_PRES_dot_AP3,P_PRES_dot_AP4,P_PRES_dot_AP5,P_PRES_dot_AP6,R_PRES_dot_P,abs_RP,breakf0; 80 double tol,Factor,sum_BREAKF,gamma,SC1,SC2,norm_of_residual,diff,L2_R,Norm_of_residual_global; 81 double *save_XPRES, *save_P_PRES, *save_R_PRES,save_R_PRES_dot_P_PRES; 82 dim_t maxit,Num_iter_global,num_iter_restart,num_iter; 83 dim_t i,z,order,n, Length_of_mem; 84 bool_t breakFlag=FALSE, maxIterFlag=FALSE, convergeFlag=FALSE,restartFlag=FALSE; 85 err_t Status=SOLVER_NO_ERROR; 86 87 88 /* adapt original routine parameters */ 89 n = Paso_SystemMatrix_getTotalNumRows(A); 90 Length_of_mem=MAX(Length_of_recursion,0)+1; 91 92 /* Test the input parameters. */ 93 if (restart>0) restart=MAX(Length_of_recursion,restart); 94 if (n < 0 || Length_of_recursion<=0) { 95 return SOLVER_INPUT_ERROR; 96 } 97 98 /* allocate memory: */ 99 X_PRES=TMPMEMALLOC(Length_of_mem,double*); 100 R_PRES=TMPMEMALLOC(Length_of_mem,double*); 101 P_PRES=TMPMEMALLOC(Length_of_mem,double*); 102 P_PRES_dot_AP=TMPMEMALLOC(Length_of_mem,double); 103 R_PRES_dot_P_PRES=TMPMEMALLOC(Length_of_mem,double); 104 BREAKF=TMPMEMALLOC(Length_of_mem,double); 105 ALPHA=TMPMEMALLOC(Length_of_mem,double); 106 AP=TMPMEMALLOC(n,double); 107 if (AP==NULL || X_PRES ==NULL || R_PRES == NULL || P_PRES == NULL || 108 P_PRES_dot_AP == NULL || R_PRES_dot_P_PRES ==NULL || BREAKF == NULL || ALPHA == NULL) { 109 Status=SOLVER_MEMORY_ERROR; 110 } else { 111 for (i=0;i6) { 280 R_PRES_dot_P=ZERO; 281 P_PRES_dot_AP0=ZERO; 282 P_PRES_dot_AP1=ZERO; 283 P_PRES_dot_AP2=ZERO; 284 P_PRES_dot_AP3=ZERO; 285 P_PRES_dot_AP4=ZERO; 286 P_PRES_dot_AP5=ZERO; 287 P_PRES_dot_AP6=ZERO; 288 #pragma omp parallel for private(z) reduction(+:R_PRES_dot_P,P_PRES_dot_AP0,P_PRES_dot_AP1,P_PRES_dot_AP2,P_PRES_dot_AP3,P_PRES_dot_AP4,P_PRES_dot_AP5,P_PRES_dot_AP6) schedule(static) 289 for (z=0;z ZERO) && (sum_BREAKF >ZERO)); 325 if (!breakFlag) { 326 breakFlag=FALSE; 327 /*** 328 ***** X_PRES and R_PRES are moved to memory: 329 ***/ 330 Factor=0.; 331 for (i=0;i0;--i) { 341 BREAKF[i]=BREAKF[i-1]; 342 R_PRES_dot_P_PRES[i]=R_PRES_dot_P_PRES[i-1]; 343 R_PRES[i]=R_PRES[i-1]; 344 X_PRES[i]=X_PRES[i-1]; 345 P_PRES[i]=P_PRES[i-1]; 346 } 347 R_PRES_dot_P_PRES[0]=save_R_PRES_dot_P_PRES; 348 R_PRES[0]=save_R_PRES; 349 X_PRES[0]=save_XPRES; 350 P_PRES[0]=save_P_PRES; 351 352 if (ABS(Factor)<=ZERO) { 353 Factor=1.; 354 BREAKF[0]=ZERO; 355 } else { 356 Factor=1./Factor; 357 BREAKF[0]=ONE; 358 } 359 for (i=0;i6) { 442 #pragma omp parallel for private(z) schedule(static) 443 for (z=0;z0.) { 468 /*** 469 ***** calculate gamma from min_(gamma){|R+gamma*(R_PRES-R)|_2}: 470 ***/ 471 SC1=ZERO; 472 SC2=ZERO; 473 L2_R=ZERO; 474 #pragma omp parallel for private(z,diff) reduction(+:SC1,SC2) schedule(static) 475 for (z=0;z0) restartFlag=(num_iter_restart >= restart); 490 } else { 491 convergeFlag=FALSE; 492 restartFlag=FALSE; 493 } 494 maxIterFlag = (num_iter >= maxit); 495 } 496 } 497 /* end of iteration */ 498 Norm_of_residual_global=norm_of_residual; 499 Num_iter_global=num_iter; 500 if (maxIterFlag) { 501 Status = SOLVER_MAXITER_REACHED; 502 } else if (breakFlag) { 503 Status = SOLVER_BREAKDOWN; 504 } 505 } 506 for (i=0;i

## Properties

Name Value
svn:eol-style native
svn:keywords Author Date Id Revision