ViewVC Help
View File | Revision Log | Show Annotations | View Changeset | Root Listing
root/freemol/trunk/smi23d/src/mengine/solve.c
Revision: 3
Committed: Mon Jun 9 21:38:26 2008 UTC (11 years, 4 months ago) by tjod
File size: 13368 byte(s)
Log Message:
test

Line File contents
1 #define EXTERN extern
2
3 // mode
4 #define NONE 0
5 #define NEWTON 1
6 #define TNCG 2
7 #define DTNCG 3
8 #define AUTO 4
9 // method
10 #define DIAG 5
11 #define SSOR 6
12 #define ICCG 7
13 #define BLOCK 8
14
15
16 #include "pcwin.h"
17 #include "pcmod.h"
18
19 #include "energies.h"
20 #include "derivs.h"
21 #include "utility.h"
22
23 //#define min(a,b) (((a) < (b) ? (a) : (b))
24
25 EXTERN struct t_minvar{
26 double cappa, stpmin, stpmax, angmax;
27 int intmax;
28 } minvar;
29
30 void solve(int,int,int,int, double *,double *,double *,double *,int *,
31 int *,int *,double *,int *,int *,int *, double (*)());
32
33
34 void precond(int,int ,int ,double *,double *,double *,int *,int *,int *,double *);
35 void cholesky(int , double *, double *);
36 void column(int,int *,int *,int *, int *,int *,int *,int *);
37
38 int stable;
39 double *f, *f_diag;
40 int *c_index, *c_value;
41
42 void solve(int mode,int method,int negtest,int nvar,double *p,double *x,
43 double *g, double *h, int *h_init, int *h_stop, int *h_index,
44 double *h_diag,int *icycle,int *iter_cg,int *fg_call,double (*fgvalue) ())
45 {
46 int i,j,k, iter,maxiter,maxhess;
47 double *m, *d, *si, *r, *q; //
48 double alpha,beta,delta,sigma;
49 double f_sigma;
50 double *x_sigma;
51 double *g_sigma;
52 double g_norm,g_rms,eps,converge;
53 double hj,gg,dq,rr,dd,rs,rs_new,r_norm;
54 int maxvar = 3*natom;
55
56 if (natom < 300)
57 maxhess = (3*natom*(3*natom-1))/2;
58 else if (natom < 800)
59 maxhess = (3*natom*(3*natom-1))/3;
60 else
61 maxhess = (3*natom*(3*natom-1))/20;
62
63
64 m = dvector(0,maxvar);
65 r = dvector(0,maxvar);
66 si = dvector(0,maxvar);
67 d = dvector(0,maxvar);
68 q = dvector(0,maxvar);
69 x_sigma = dvector(0,maxvar);
70 g_sigma = dvector(0,maxvar);
71 f_diag = dvector(0,maxvar);
72 f = dvector(0,maxhess);
73 c_index = ivector(0,maxhess);
74 c_value = ivector(0,maxhess);
75
76 if (mode != DTNCG && method != NONE)
77 {
78 for (i=0; i < nvar; i++)
79 m[i] = 1.0 / sqrt(fabs(h_diag[i]));
80
81 for (i=0; i < nvar; i++)
82 {
83 g[i] = g[i] * m[i];
84 h_diag[i] *= m[i] * m[i];
85 for (j = h_init[i]; j < h_stop[i]; j++)
86 {
87 k = h_index[j];
88 h[j] *= m[i] * m[k];
89 }
90 }
91 }
92
93 iter = 0;
94 gg = 0.0;
95 for (i=0; i < nvar; i++)
96 {
97 p[i] = 0.0;
98 r[i] = -g[i];
99 gg += g[i]*g[i];
100 }
101 g_norm = sqrt(gg);
102 precond(method, iter,nvar, si, r, h, h_init, h_stop, h_index, h_diag);
103 rs = 0.0;
104 for (i=0; i < nvar; i++)
105 {
106 d[i] = si[i];
107 rs += r[i]*si[i];
108 }
109
110 if (mode == NEWTON)
111 {
112 eps = 1.0e-10;
113 maxiter = nvar;
114 } else if (mode == TNCG || mode == DTNCG)
115 {
116 delta = 1.0;
117 eps = delta/(double) *icycle;
118 g_rms = g_norm/sqrt((double) nvar);
119 eps = MinFun(eps,g_rms);
120 converge = 1.0;
121 maxiter = (int) (10.0*sqrt((double)nvar));
122 }
123 iter = 1;
124
125 // evaluate the matrix vector product
126 while (TRUE)
127 {
128 if (mode == TNCG || mode == NEWTON)
129 {
130 for (i=0; i < nvar; i++)
131 q[i] = 0.0;
132 for (i=0; i < nvar; i++)
133 {
134 q[i] += h_diag[i]*d[i];
135 for(j= h_init[i]; j < h_stop[i]; j++)
136 {
137 k = h_index[j];
138 hj = h[j];
139 q[i] += hj*d[k];
140 q[k] += hj*d[i];
141 }
142 }
143 }else if (mode == DTNCG)
144 {
145 dd = 0.0;
146 for (i=0; i < nvar; i++)
147 dd += d[i]*d[i];
148
149 sigma = 1.0e-7 / sqrt(dd);
150
151 for (i=0; i < nvar; i++)
152 x_sigma[i] = x[i] + sigma*d[i];
153
154 fg_call = fg_call + 1;
155 f_sigma = fgvalue (x_sigma,g_sigma);
156
157 for (i=0; i < nvar; i++)
158 q[i] = (g_sigma[i]-g[i]) / sigma;
159 }
160
161 // test for direction of negative curvature
162 dq = 0.0;
163 for (i=0; i < nvar; i++)
164 dq += d[i]*q[i];
165
166 if (negtest)
167 {
168 if (dq <= 0.0)
169 {
170 if (iter == 1)
171 {
172 for (i=0; i < nvar; i++)
173 p[i] = d[i];
174 }
175 goto L_10;
176 }
177 }
178
179 // test truncated Newton termination criterion
180 alpha = rs / dq;
181 rr = 0.0;
182 for (i=0; i < nvar; i++)
183 {
184 p[i] += alpha*d[i];
185 r[i] -= alpha*q[i];
186 rr += r[i]*r[i];
187 }
188 r_norm = sqrt(rr);
189 if (r_norm/g_norm <= eps)
190 goto L_10;
191
192 // precondition
193 precond(method,iter,nvar, si, r,h,h_init,h_stop,
194 h_index,h_diag);
195
196 // update direction
197 rs_new = 0.0;
198 for (i=0; i < nvar; i++)
199 rs_new += r[i]*si[i];
200
201 beta = rs_new/rs;
202 rs = rs_new;
203
204 for (i=0; i < nvar; i++)
205 d[i] = si[i] + beta*d[i];
206
207 // check for overlimit next iteration
208 if (iter >= maxiter)
209 goto L_10;
210
211 iter++;
212 }
213
214 // termination
215
216 L_10:
217 if (mode != DTNCG && method != NONE)
218 {
219 for (i=0; i < nvar; i++)
220 {
221 p[i] *= m[i];
222 g[i] /= m[i];
223 }
224 }
225 iter_cg += iter;
226 free_dvector( m ,0,maxvar);
227 free_dvector( r ,0,maxvar);
228 free_dvector( si ,0,maxvar);
229 free_dvector( d ,0,maxvar);
230 free_dvector( q ,0,maxvar);
231 free_dvector( x_sigma ,0,maxvar);
232 free_dvector( g_sigma ,0,maxvar);
233 free_dvector( f_diag ,0,maxvar);
234 free_dvector( f ,0,maxhess);
235 free_ivector( c_index ,0,maxhess);
236 free_ivector( c_value ,0,maxhess);
237 }
238
239 void precond(int method, int iter, int nvar, double *s, double *r,
240 double *h, int *h_init, int *h_stop,int *h_index,
241 double *h_diag)
242 {
243 int i,j,k,ii,kk,iii,kkk,nblock,ix,iy,iz,icount;
244 int *c_init, *c_stop; // c_init[maxvar],c_stop[maxvar];
245 double a[6],b[3];
246 double omega,factor,*diag; // diag[maxvar];
247 double maxalpha, alpha, f_i, f_k;
248 int maxvar;
249
250 maxvar = 3*natom;
251 c_init = ivector(0,maxvar);
252 c_stop = ivector(0,maxvar);
253 diag = dvector(0,maxvar);
254
255
256 // no preconditioning
257 if (method == NONE)
258 {
259 for (i=0; i < nvar; i++)
260 s[i] = r[i];
261 }
262 // diagonal
263 if (method == DIAG)
264 {
265 for (i=0; i < nvar; i++)
266 s[i] = r[i]/ fabs(h_diag[i]);
267 }
268 // block diagonal
269 if (method == BLOCK)
270 {
271 nblock = 3;
272 for (i=1; i <= natom; i++)
273 {
274 iz = (i-1)*3+2;
275 iy = (i-1)*3+1;
276 ix = (i-1)*3;
277 a[0] = h_diag[ix];
278 if (h_index[h_init[ix]] == iy)
279 a[1] = h[h_init[ix]];
280 else
281 a[1] = 0.0;
282
283 if (h_index[h_init[ix]+1] == iz)
284 a[2] = h[h_init[ix]+1];
285 else
286 a[2] = 0.0;
287
288 a[3] = h_diag[iy];
289 if (h_index[h_init[iy]] == iz)
290 a[4] = h[h_init[iy]];
291 else
292 a[4] = 0.0;
293
294 a[5] = h_diag[iz];
295 b[0] = r[ix];
296 b[1] = r[iy];
297 b[2] = r[iz];
298 cholesky (nblock,a,b);
299 s[ix] = b[0];
300 s[iy] = b[1];
301 s[iz] = b[2];
302 }
303 }
304
305 // SSOR
306 if (method == SSOR)
307 {
308 omega = 1.0;
309 factor = 2.0 - omega;
310 for (i=0; i < nvar; i++)
311 {
312 s[i] = r[i] * factor;
313 diag[i] = h_diag[i] / omega;
314 }
315 for (i=0; i < nvar; i++)
316 {
317 s[i] = s[i] / diag[i];
318 for (j=h_init[i]; j< h_stop[i]; j++)
319 {
320 k = h_index[j];
321 s[k] = s[k] - h[j]*s[i];
322 }
323 }
324 for (i=nvar-1; i >= 0; i--)
325 {
326 s[i] = s[i] * diag[i];
327 for (j=h_init[i]; j< h_stop[i]; j++)
328 {
329 k = h_index[j];
330 s[i] = s[i] - h[j]*s[k];
331 }
332 s[i] = s[i] / diag[i];
333 }
334 }
335 // incomplete cholesky preconditioning
336 if (method == ICCG && iter == 0)
337 {
338 column (nvar,h_init,h_stop,h_index,
339 c_init,c_stop,c_index,c_value);
340 stable = TRUE;
341 icount = 0;
342 maxalpha = 2.1;
343 alpha = -0.001;
344 L_10:
345 if (alpha <= 0.0)
346 alpha = alpha + 0.001;
347 else
348 alpha = 2.0 * alpha;
349
350 if (alpha > maxalpha)
351 {
352 stable = FALSE;
353 fprintf(pcmoutfile,"Precond - Incomplete Cholesky is Unstable\n");
354 } else
355 {
356 factor = 1.0 + alpha;
357 for (i=0; i < nvar; i++)
358 {
359 f_diag[i] = factor * h_diag[i];
360 for(j = c_init[i]; j<= c_stop[i]; j++)
361 {
362 k = c_index[j];
363 f_i = f[c_value[j]];
364 f_diag[i] = f_diag[i] - (f_i*f_i*f_diag[k]);
365 icount++;
366 }
367 if (f_diag[i] <= 0.0) goto L_10;
368 if (f_diag[i] < 1.0e-7) f_diag[i] = 1.0e-7;
369 f_diag[i] = 1.0 / f_diag[i];
370 for( j = h_init[i]; j < h_stop[i]; j++)
371 {
372 k = h_index[j];
373 f[j] = h[j];
374 ii = c_init[i];
375 kk = c_init[k];
376 while (ii <= c_stop[i] && kk <= c_stop[k])
377 {
378 iii = c_index[ii];
379 kkk = c_index[kk];
380 if (iii < kkk)
381 ii = ii + 1;
382 else if (kkk < iii)
383 kk = kk + 1;
384 else
385 {
386 f_i = f[c_value[ii]];
387 f_k = f[c_value[kk]];
388 f[j] = f[j] - (f_i*f_k*f_diag[iii]);
389 ii = ii + 1;
390 kk = kk + 1;
391 icount = icount + 1;
392 }
393 }
394 }
395 }
396 }
397 }
398 // incomplete cholesky preconditioning (solution phase)
399
400 if (method == ICCG)
401 {
402 if (stable)
403 {
404 for (i=0; i < nvar; i++)
405 s[i] = r[i];
406
407 for (i=0; i < nvar; i++)
408 {
409 s[i] = s[i] * f_diag[i];
410 for( j = h_init[i]; j < h_stop[i]; j++)
411 {
412 k = h_index[j];
413 s[k] = s[k] - f[j]*s[i];
414 }
415 }
416 for (i = nvar-1; i >= 0; i--)
417 {
418 s[i] = s[i] / f_diag[i];
419 for( j = h_init[i]; j < h_stop[i]; j++)
420 {
421 k = h_index[j];
422 s[i] = s[i] - f[j]*s[k];
423 }
424 s[i] = s[i] * f_diag[i];
425 }
426 } else
427 {
428 for (i=0; i < nvar; i++)
429 s[i] = r[i] / fabs(h_diag[i]);
430 }
431 }
432 free_ivector(c_init ,0,maxvar);
433 free_ivector(c_stop ,0,maxvar);
434 free_dvector(diag ,0,maxvar);
435 }
436
437 void cholesky(int nvar, double *a, double *b)
438 {
439 int i,j,k,ii,ij,ik,im,jk,jm,ki,kk;
440 double r,sa,t;
441
442 ii = 0;
443 for (i=0; i < nvar; i++)
444 {
445 im = i - 1;
446 if (i != 0)
447 {
448 ij = i;
449 for( j = 0; j < im; j++)
450 {
451 r = a[ij];
452 if (j != 0)
453 {
454 ik = i;
455 jk = j;
456 jm = j - 1;
457 for ( k = 0; k < jm; k++)
458 {
459 r = r - a[ik]*a[jk];
460 ik = nvar-1 - k + ik;
461 jk = nvar-1 - k + jk;
462 }
463 }
464 a[ij] = r;
465 ij = nvar-1 - j + ij;
466 }
467 }
468 r = a[ii];
469 if (i != 0)
470 {
471 kk = 0;
472 ik = i;
473 for (k = 0; k < im; k++)
474 {
475 sa = a[ik];
476 t = sa * a[kk];
477 a[ik] = t;
478 r = r - sa*t;
479 ik = nvar-1 - k + ik;
480 kk = nvar-1 - k +1 + kk;
481 }
482 }
483 a[ii] = 1.0/ r;
484 ii = nvar-1 - i + 1 + ii;
485 }
486
487 // solve linear equations; first solve Ly = b for y
488
489 for (i=0; i < nvar; i++)
490 {
491 if (i != 0)
492 {
493 ik = i;
494 im = i - 1;
495 r = b[i];
496 for (k = 0; k < im; k++)
497 {
498 r = r - b[k]*a[ik];
499 ik = nvar-1 - k + ik;
500 }
501 b[i] = r;
502 }
503 }
504
505 // now, solve (D)(L transpose)(x) = y for x
506
507 ii = (nvar*(nvar+1)/2)-1;
508 for ( j = 0; j < nvar; j++)
509 {
510 i = nvar-1 - j;
511 r = b[i] * a[ii];
512 if (j != 0)
513 {
514 im = i + 1;
515 ki = ii + 1;
516 for( k = im; k < nvar; k++)
517 {
518 r = r - a[ki]*b[k];
519 ki++;
520 }
521 }
522 b[i] = r;
523 ii = ii - j - 2;
524 }
525 }