Annotation of /trunk/paso/src/GMRES2.c

Revision 1798 - (hide annotations)
Wed Sep 17 06:21:12 2008 UTC (13 years ago) by gross
File MIME type: text/plain
File size: 6604 byte(s)
```Fixes for the JacobeanFreeNewton scheme. Still needs to be tested under OPENMP but runs under MPI.

```
 1 gross 1476 /* \$Id:\$ */ 2 3 /******************************************************* 4 * 5 * Copyright 2008 by University of Queensland 6 * 7 * http://esscc.uq.edu.au 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 * Purpose 15 * ======= 16 * 17 gross 1639 * GMRES2 solves the non-linear system f0+J_0*d=0 18 gross 1476 * where f0=F(x0), J_0 is the jacobian of F at x0. 19 * 20 * Convergence test: norm(f0+J_0*d)<=tolerance*norm(f0) 21 * 22 * Arguments 23 * ========= 24 * 25 * Paso_Function * F evaluation of F (including any preconditioner) 26 * 27 * x0 (input) current point 28 * 29 * f0 (input) function F at current point x0 30 * 31 * d (output) solution of f0+J0*d=0 with accuracy tolerance 32 * 33 * iter (input/output) 34 * On input, the maximum num_iterations to be performed. 35 * On output, actual number of num_iterations performed. 36 * tolerance (input/output) 37 * On input, the allowable convergence measure for norm(f0+J0*d)/norm(f0) 38 * On output, the final value for norm(f0+J0*d) 39 * return value 40 * 41 * =SOLVER_NO_ERROR: Successful exit. approximate solution returned. 42 * =SOLVER_MAXNUM_ITER_REACHED 43 * =SOLVER_INPUT_ERROR Illegal parameter: 44 * =SOLVER_BREAKDOWN: bad luck! 45 * =SOLVER_MEMORY_ERROR : no memory available 46 * 47 * ============================================================== 48 */ 49 50 #include "Common.h" 51 #include "Solver.h" 52 gross 1639 #include "PasoUtil.h" 53 gross 1476 #ifdef _OPENMP 54 #include 55 #endif 56 57 gross 1639 err_t Paso_Solver_GMRES2( 58 gross 1476 Paso_Function * F, 59 const double* f0, 60 const double* x0, 61 double * x, 62 dim_t *iter, 63 double* tolerance, 64 Paso_Performance* pp) 65 { 66 double static RENORMALIZATION_CONST=0.001; 67 gross 1639 dim_t l=(*iter)+1, iter_max=*iter,k=0,i,j; 68 const dim_t n=F->n; 69 gross 1476 double rel_tol=*tolerance, abs_tol, normf0, normv, normv2, hh, hr, nu, norm_of_residual=0.; 70 bool_t breakFlag=FALSE, maxIterFlag=FALSE, convergeFlag=FALSE; 71 double *h=NULL, **v=NULL, *c=NULL,*s=NULL,*g=NULL, *work=NULL; 72 err_t Status=SOLVER_NO_ERROR; 73 74 /* Test the input parameters. */ 75 if (n < 0 || iter_max<=0 || l<1 || rel_tol<0) { 76 return SOLVER_INPUT_ERROR; 77 } 78 Paso_zeroes(n,x); 79 /* 80 * allocate memory: 81 */ 82 h=TMPMEMALLOC(l*l,double); 83 v=TMPMEMALLOC(iter_max,double*); 84 c=TMPMEMALLOC(l,double); 85 s=TMPMEMALLOC(l,double); 86 g=TMPMEMALLOC(l,double); 87 work=TMPMEMALLOC(n,double); 88 89 if (h==NULL || v ==NULL || c== NULL || s == NULL || g== NULL || work==NULL) { 90 Status=SOLVER_MEMORY_ERROR; 91 } 92 if ( Status ==SOLVER_NO_ERROR ) { 93 gross 1639 for (i=0;impi_info); 98 k=0; 99 convergeFlag=(ABS(normf0)<=0); 100 if (! convergeFlag) { 101 abs_tol=rel_tol*normf0; 102 gross 1639 v[0]=TMPMEMALLOC(n,double); 103 if (v[0]==NULL) { 104 Status=SOLVER_MEMORY_ERROR; 105 } else { 106 Paso_zeroes(n,v[0]); 107 Paso_Update(n,1.,v[0],-1./normf0,f0); 108 g[0]=normf0; 109 } 110 gross 1798 while( (! (breakFlag || maxIterFlag || convergeFlag)) && (Status ==SOLVER_NO_ERROR) ) { 111 gross 1476 k++; 112 gross 1639 v[k]=TMPMEMALLOC(n,double); 113 if (v[k]==NULL) { 114 Status=SOLVER_MEMORY_ERROR; 115 } else { 116 /* 117 * call directional derivative function 118 */ 119 gross 1798 Paso_FunctionDerivative(v[k],v[k-1],F,f0,x0,work,TRUE); 120 gross 1639 normv=Paso_l2(n,v[k],F->mpi_info); 121 /* 122 * Modified Gram-Schmidt 123 */ 124 for (j=0;jmpi_info); 126 Paso_Update(n,1.,v[k],(-hh),v[j]); 127 h[INDEX2(j,k-1,l)]=hh; 128 } 129 normv2=Paso_l2(n,v[k],F->mpi_info); 130 h[INDEX2(k,k-1,l)]=normv2; 131 /* 132 * reorthogonalize 133 */ 134 if (! (normv + RENORMALIZATION_CONST*normv2 > normv)) { 135 for (j=0;jmpi_info); 137 h[INDEX2(j,k-1,l)]+=hr; 138 Paso_Update(n,1.,v[k],(-hr),v[j]); 139 } 140 normv2=Paso_l2(n,v[k],F->mpi_info); 141 h[INDEX2(k,k-1,l)]=normv2; 142 } 143 /* 144 * watch out for happy breakdown 145 */ 146 if(normv2 > 0.) { 147 Paso_Update(n,1./normv2,v[k],0.,v[k]); 148 } 149 /* 150 * Form and store the information for the new Givens rotation 151 */ 152 ApplyGivensRotations(k,&h[INDEX2(0,k-1,l)],c,s); 153 /* 154 * Don't divide by zero if solution has been found 155 */ 156 g[k]=0; 157 nu=sqrt(h[INDEX2(k-1,k-1,l)]*h[INDEX2(k-1,k-1,l)]+h[INDEX2(k,k-1,l)]*h[INDEX2(k,k-1,l)]); 158 if (nu>0) { 159 c[k-1]=h[INDEX2(k-1,k-1,l)]/nu; 160 gross 1798 s[k-1]=-h[INDEX2(k,k-1,l)]/nu; 161 h[INDEX2(k-1,k-1,l)]=c[k-1]*h[INDEX2(k-1,k-1,l)]-s[k-1]*h[INDEX2(k,k-1,l)]; 162 gross 1639 h[INDEX2(k,k-1,l)]=0; 163 ApplyGivensRotations(2,&(g[k-1]),&(c[k-1]),&(s[k-1])); 164 } 165 norm_of_residual=fabs(g[k]); 166 maxIterFlag = (k>=iter_max); 167 convergeFlag = (fabs(g[k]) <= abs_tol); 168 printf("GMRES2 step %d: error %e (tol=%e)\n",k,fabs(g[k]),abs_tol); 169 } 170 gross 1476 } 171 } 172 /* 173 * all done and ready for the forward substitution: 174 */ 175 gross 1798 for (i=k-1;i>=0;--i) { 176 for (j=i+1;j