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

Revision 2108 - (show annotations)
Fri Nov 28 05:09:23 2008 UTC (13 years, 7 months ago) by gross
File MIME type: text/plain
File size: 27271 byte(s)
```some minor changes to PCG and some extra suspicious characters.
```
 1 2 /******************************************************* 3 * 4 * Copyright (c) 2003-2008 by University of Queensland 5 * Earth Systems Science Computational Center (ESSCC) 6 * http://www.uq.edu.au/esscc 7 * 8 * Primary Business: Queensland, Australia 9 * Licensed under the Open Software License version 3.0 10 * http://www.opensource.org/licenses/osl-3.0.php 11 * 12 *******************************************************/ 13 14 15 /* 16 * Purpose 17 * ======= 18 * 19 * GMRES solves the linear system A*x=b using the 20 * truncated and restered GMRES method with preconditioning. 21 * 22 * Convergence test: norm( b - A*x )< TOL. 23 * 24 * 25 * 26 * Arguments 27 * ========= 28 * 29 * r (input/output) double array, dimension n. 30 * On entry, residual of inital guess X 31 * 32 * x (input/output) double array, dimension n. 33 * On input, the initial guess. 34 * 35 * iter (input/output) int 36 * On input, the maximum num_iterations to be performed. 37 * On output, actual number of num_iterations performed. 38 * 39 * tolerance (input/output) DOUBLE PRECISION 40 * On input, the allowable convergence measure for 41 * norm( b - A*x ) 42 * On output, the final value of this measure. 43 * 44 * Length_of_recursion (input) gives the number of residual to be kept in orthogonalization process 45 * 46 * restart (input) If restart>0, iteration is resterted a after restart steps. 47 * 48 * INFO (output) int 49 * 50 * =SOLVER_NO_ERROR: Successful exit. num_iterated approximate solution returned. 51 * =SOLVER_MAXNUM_ITER_REACHED 52 * =SOLVER_INPUT_ERROR Illegal parameter: 53 * =SOLVER_BREAKDOWN: bad luck! 54 * =SOLVER_MEMORY_ERROR : no memory available 55 * 56 * ============================================================== 57 */ 58 59 #include "Common.h" 60 #include "SystemMatrix.h" 61 #include "Solver.h" 62 #ifdef _OPENMP 63 #include 64 #endif 65 66 err_t Paso_Solver_GMRES( 67 Paso_SystemMatrix * A, 68 double * r, 69 double * x, 70 dim_t *iter, 71 double * tolerance,dim_t Length_of_recursion,dim_t restart, 72 Paso_Performance* pp) { 73 74 /* Local variables */ 75 76 #ifdef _OPENMP 77 const int num_threads=omp_get_max_threads(); 78 #else 79 const int num_threads=1; 80 #endif 81 double *AP,**X_PRES,**R_PRES,**P_PRES, *dots, *loc_dots; 82 double *P_PRES_dot_AP,*R_PRES_dot_P_PRES,*BREAKF,*ALPHA; 83 double R_PRES_dot_AP0,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,breakf0; 84 double tol,Factor,sum_BREAKF,gamma,SC1,SC2,norm_of_residual=0,diff,L2_R,Norm_of_residual_global=0; 85 double *save_XPRES, *save_P_PRES, *save_R_PRES,save_R_PRES_dot_P_PRES; 86 dim_t maxit,Num_iter_global=0,num_iter_restart=0,num_iter; 87 dim_t i,z,order,n, Length_of_mem, th, local_n , rest, n_start ,n_end; 88 bool_t breakFlag=FALSE, maxIterFlag=FALSE, convergeFlag=FALSE,restartFlag=FALSE; 89 err_t Status=SOLVER_NO_ERROR; 90 91 92 /* adapt original routine parameters */ 93 n = Paso_SystemMatrix_getTotalNumRows(A); 94 Length_of_mem=MAX(Length_of_recursion,0)+1; 95 96 /* Test the input parameters. */ 97 if (restart>0) restart=MAX(Length_of_recursion,restart); 98 if (n < 0 || Length_of_recursion<=0) { 99 return SOLVER_INPUT_ERROR; 100 } 101 102 /* allocate memory: */ 103 X_PRES=TMPMEMALLOC(Length_of_mem,double*); 104 R_PRES=TMPMEMALLOC(Length_of_mem,double*); 105 P_PRES=TMPMEMALLOC(Length_of_mem,double*); 106 loc_dots=TMPMEMALLOC(MAX(Length_of_mem+1,3),double); 107 dots=TMPMEMALLOC(MAX(Length_of_mem+1,3),double); 108 P_PRES_dot_AP=TMPMEMALLOC(Length_of_mem,double); 109 R_PRES_dot_P_PRES=TMPMEMALLOC(Length_of_mem,double); 110 BREAKF=TMPMEMALLOC(Length_of_mem,double); 111 ALPHA=TMPMEMALLOC(Length_of_mem,double); 112 AP=TMPMEMALLOC(n,double); 113 if (AP==NULL || X_PRES ==NULL || R_PRES == NULL || P_PRES == NULL || 114 P_PRES_dot_AP == NULL || R_PRES_dot_P_PRES ==NULL || BREAKF == NULL || ALPHA == NULL || dots==NULL || loc_dots==NULL) { 115 Status=SOLVER_MEMORY_ERROR; 116 } else { 117 for (i=0;impi_info->comm); 366 R_PRES_dot_P_PRES[0]=dots[0]; 367 memcpy(P_PRES_dot_AP,&dots[1],sizeof(double)*order); 368 #else 369 R_PRES_dot_P_PRES[0]=loc_dots[0]; 370 memcpy(P_PRES_dot_AP,&loc_dots[1],sizeof(double)*order); 371 #endif 372 R_PRES_dot_AP0=R_PRES_dot_P_PRES[0]; 373 374 /*** if sum_BREAKF is equal to zero a breakdown occurs. 375 *** iteration procedure can be continued but R_PRES is not the 376 *** residual of X_PRES approximation. 377 ***/ 378 sum_BREAKF=0.; 379 #pragma ivdep 380 for (i=0;i PASO_ZERO) && (sum_BREAKF >PASO_ZERO)); 382 if (!breakFlag) { 383 breakFlag=FALSE; 384 /*** 385 ***** X_PRES and R_PRES are moved to memory: 386 ***/ 387 Factor=0.; 388 #pragma ivdep 389 for (i=0;i0;--i) { 400 BREAKF[i]=BREAKF[i-1]; 401 R_PRES_dot_P_PRES[i]=R_PRES_dot_P_PRES[i-1]; 402 R_PRES[i]=R_PRES[i-1]; 403 X_PRES[i]=X_PRES[i-1]; 404 P_PRES[i]=P_PRES[i-1]; 405 } 406 R_PRES_dot_P_PRES[0]=save_R_PRES_dot_P_PRES; 407 R_PRES[0]=save_R_PRES; 408 X_PRES[0]=save_XPRES; 409 P_PRES[0]=save_P_PRES; 410 411 if (ABS(Factor)<=PASO_ZERO) { 412 Factor=1.; 413 BREAKF[0]=PASO_ZERO; 414 } else { 415 Factor=1./Factor; 416 BREAKF[0]=PASO_ONE; 417 } 418 for (i=0;i0.) { 534 /*** 535 ***** calculate gamma from min_(gamma){|R+gamma*(R_PRES-R)|_2}: 536 ***/ 537 memset(loc_dots,0,sizeof(double)*3); 538 #pragma omp parallel for private(th,z,i,local_n, rest, n_start, n_end,diff,SC1,SC2) 539 for (th=0;thmpi_info->comm); 560 SC1=dots[0]; 561 SC2=dots[1]; 562 #else 563 SC1=loc_dots[0]; 564 SC2=loc_dots[1]; 565 #endif 566 gamma=(SC1<=PASO_ZERO) ? PASO_ZERO : -SC2/SC1; 567 #pragma omp parallel for private(th,z,local_n, rest, n_start, n_end,diff,L2_R) 568 for (th=0;thmpi_info->comm); 587 L2_R=dots[2]; 588 #else 589 L2_R=loc_dots[2]; 590 #endif 591 norm_of_residual=sqrt(L2_R); 592 convergeFlag = (norm_of_residual <= tol); 593 if (restart>0) restartFlag=(num_iter_restart >= restart); 594 } else { 595 convergeFlag=FALSE; 596 restartFlag=FALSE; 597 } 598 maxIterFlag = (num_iter >= maxit); 599 } 600 } 601 /* end of iteration */ 602 Norm_of_residual_global=norm_of_residual; 603 Num_iter_global=num_iter; 604 if (maxIterFlag) { 605 Status = SOLVER_MAXITER_REACHED; 606 } else if (breakFlag) { 607 Status = SOLVER_BREAKDOWN; 608 } 609 } 610 for (i=0;i

## Properties

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