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

Revision 2548 - (show annotations)
Mon Jul 20 06:20:06 2009 UTC (10 years, 3 months ago) by jfenwick
File MIME type: text/plain
File size: 27270 byte(s)
```Updating copyright notices
```
 1 2 /******************************************************* 3 * 4 * Copyright (c) 2003-2009 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 /*** if sum_BREAKF is equal to zero a breakdown occurs. 374 *** iteration procedure can be continued but R_PRES is not the 375 *** residual of X_PRES approximation. 376 ***/ 377 sum_BREAKF=0.; 378 #pragma ivdep 379 for (i=0;i PASO_ZERO) && (sum_BREAKF >PASO_ZERO)); 381 if (!breakFlag) { 382 breakFlag=FALSE; 383 /*** 384 ***** X_PRES and R_PRES are moved to memory: 385 ***/ 386 Factor=0.; 387 #pragma ivdep 388 for (i=0;i0;--i) { 399 BREAKF[i]=BREAKF[i-1]; 400 R_PRES_dot_P_PRES[i]=R_PRES_dot_P_PRES[i-1]; 401 R_PRES[i]=R_PRES[i-1]; 402 X_PRES[i]=X_PRES[i-1]; 403 P_PRES[i]=P_PRES[i-1]; 404 } 405 R_PRES_dot_P_PRES[0]=save_R_PRES_dot_P_PRES; 406 R_PRES[0]=save_R_PRES; 407 X_PRES[0]=save_XPRES; 408 P_PRES[0]=save_P_PRES; 409 410 if (ABS(Factor)<=PASO_ZERO) { 411 Factor=1.; 412 BREAKF[0]=PASO_ZERO; 413 } else { 414 Factor=1./Factor; 415 BREAKF[0]=PASO_ONE; 416 } 417 for (i=0;i0.) { 533 /*** 534 ***** calculate gamma from min_(gamma){|R+gamma*(R_PRES-R)|_2}: 535 ***/ 536 memset(loc_dots,0,sizeof(double)*3); 537 #pragma omp parallel for private(th,z,i,local_n, rest, n_start, n_end,diff,SC1,SC2) 538 for (th=0;thmpi_info->comm); 559 SC1=dots[0]; 560 SC2=dots[1]; 561 #else 562 SC1=loc_dots[0]; 563 SC2=loc_dots[1]; 564 #endif 565 gamma=(SC1<=PASO_ZERO) ? PASO_ZERO : -SC2/SC1; 566 #pragma omp parallel for private(th,z,local_n, rest, n_start, n_end,diff,L2_R) 567 for (th=0;thmpi_info->comm); 586 L2_R=dots[2]; 587 #else 588 L2_R=loc_dots[2]; 589 #endif 590 norm_of_residual=sqrt(L2_R); 591 convergeFlag = (norm_of_residual <= tol); 592 if (restart>0) restartFlag=(num_iter_restart >= restart); 593 } else { 594 convergeFlag=FALSE; 595 restartFlag=FALSE; 596 } 597 maxIterFlag = (num_iter >= maxit); 598 } 599 } 600 /* end of iteration */ 601 Norm_of_residual_global=norm_of_residual; 602 Num_iter_global=num_iter; 603 if (maxIterFlag) { 604 Status = SOLVER_MAXITER_REACHED; 605 } else if (breakFlag) { 606 Status = SOLVER_BREAKDOWN; 607 } 608 } 609 for (i=0;i

## Properties

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