# Contents of /branches/doubleplusgood/paso/src/GMRES.cpp

Revision 4291 - (show annotations)
Thu Mar 7 08:57:58 2013 UTC (6 years, 1 month ago) by jfenwick
File size: 27247 byte(s)

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

## Properties

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