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

Revision 1388 - (show annotations)
Fri Jan 11 07:45:58 2008 UTC (11 years, 8 months ago) by trankine
File MIME type: text/plain
File size: 24161 byte(s)
```And get the *(&(*&(* name right
```
 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) { 331 #pragma omp master 332 { 333 R_PRES_dot_P=ZERO; 334 P_PRES_dot_AP0=ZERO; 335 P_PRES_dot_AP1=ZERO; 336 P_PRES_dot_AP2=ZERO; 337 P_PRES_dot_AP3=ZERO; 338 P_PRES_dot_AP4=ZERO; 339 P_PRES_dot_AP5=ZERO; 340 P_PRES_dot_AP6=ZERO; 341 } 342 #pragma omp barrier 343 #pragma omp 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) 344 for (z=0;z ZERO) && (sum_BREAKF >ZERO)); 389 if (!breakFlag) { 390 breakFlag=FALSE; 391 /*** 392 ***** X_PRES and R_PRES are moved to memory: 393 ***/ 394 #pragma omp master 395 { 396 Factor=0.; 397 for (i=0;i0;--i) { 407 BREAKF[i]=BREAKF[i-1]; 408 R_PRES_dot_P_PRES[i]=R_PRES_dot_P_PRES[i-1]; 409 R_PRES[i]=R_PRES[i-1]; 410 X_PRES[i]=X_PRES[i-1]; 411 P_PRES[i]=P_PRES[i-1]; 412 } 413 R_PRES_dot_P_PRES[0]=save_R_PRES_dot_P_PRES; 414 R_PRES[0]=save_R_PRES; 415 X_PRES[0]=save_XPRES; 416 P_PRES[0]=save_P_PRES; 417 418 if (ABS(Factor)<=ZERO) { 419 Factor=1.; 420 BREAKF[0]=ZERO; 421 } else { 422 Factor=1./Factor; 423 BREAKF[0]=ONE; 424 } 425 for (i=0;i6) { 510 #pragma omp for private(z) schedule(static) 511 for (z=0;z0.) { 536 /*** 537 ***** calculate gamma from min_(gamma){|R+gamma*(R_PRES-R)|_2}: 538 ***/ 539 #pragma omp master 540 { 541 SC1=ZERO; 542 SC2=ZERO; 543 L2_R=ZERO; 544 } 545 #pragma omp barrier 546 #pragma omp for private(z,diff) reduction(+:SC1,SC2) schedule(static) 547 for (z=0;z0) restartFlag=(num_iter_restart >= restart); 562 } else { 563 convergeFlag=FALSE; 564 restartFlag=FALSE; 565 } 566 maxIterFlag = (num_iter >= maxit); 567 } 568 } 569 /* end of iteration */ 570 #pragma omp master 571 { 572 Norm_of_residual_global=norm_of_residual; 573 Num_iter_global=num_iter; 574 if (maxIterFlag) { 575 Status = SOLVER_MAXITER_REACHED; 576 } else if (breakFlag) { 577 Status = SOLVER_BREAKDOWN; 578 } 579 } 580 } /* end of parallel region */ 581 for (i=0;i

## Properties

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