 |
Millepede-II
V04-00-00
|
|
Millepede-II
V04-00-00
|
00001 00002 ! Code converted using TO_F90 by Alan Miller 00003 ! Date: 2012-03-16 Time: 11:06:43 00004 00007 00288 00289 SUBROUTINE minres( n, b, r1, r2, v, w, w1, w2, x, y, & 00290 aprod, msolve, checka, precon, shift, nout , itnlim, rtol, & 00291 istop, itn, anorm, acond, rnorm, ynorm ) 00292 00293 IMPLICIT NONE 00294 EXTERNAL aprod, msolve 00295 INTEGER :: n, nout, itnlim, istop, itn 00296 LOGICAL :: checka, precon 00297 DOUBLE PRECISION :: shift, rtol, anorm, acond, rnorm, ynorm, 00298 b(n), r1(n), r2(n), v(n), w(n), w1(n), w2(n), x(n), y(n) 00299 00300 EXTERNAL ddot , dnrm2 00301 DOUBLE PRECISION :: ddot , dnrm2 00302 00303 ! Local variables 00304 00305 DOUBLE PRECISION :: alfa , beta , beta1 , cs , 00306 dbar , delta , denom , diag , eps , epsa , epsln , epsr , epsx , 00307 agamma, gbar , gmax , gmin , oldb , oldeps, qrnorm, phi , phibar, 00308 rhs1 , rhs2 , s , sn , t , tnorm2, ynorm2, z 00309 INTEGER :: i 00310 LOGICAL :: debug, prnt 00311 00312 DOUBLE PRECISION :: zero, one, two, ten 00313 PARAMETER ( zero = 0.0D+0, one = 1.0D+0, & 00314 two = 2.0D+0, ten = 10.0D+0 ) 00315 00316 CHARACTER (LEN=16) :: enter, EXIT 00317 CHARACTER (LEN=52) :: msg(-1:8) 00318 00319 DATA enter /' Enter MINRES. '/, EXIT /' Exit MINRES. '/ 00320 00321 DATA msg & 00322 / 'beta2 = 0. If M = I, b and x are eigenvectors of A', & 00323 'beta1 = 0. The exact solution is x = 0', & 00324 'Requested accuracy achieved, as determined by rtol', & 00325 'Reasonable accuracy achieved, given eps', & 00326 'x has converged to an eigenvector', 'Acond has exceeded 0.1/eps', & 00327 'The iteration limit was reached', & 00328 'Aprod does not define a symmetric matrix', & 00329 'Msolve does not define a symmetric matrix', & 00330 'Msolve does not define a pos-def preconditioner' / 00331 ! ------------------------------------------------------------------ 00332 00333 debug = .false. 00334 00335 ! ------------------------------------------------------------------ 00336 ! Compute eps, the machine precision. The call to daxpy is 00337 ! intended to fool compilers that use extra-length registers. 00338 ! 31 May 1999: Hardwire eps so the debugger can step thru easily. 00339 ! ------------------------------------------------------------------ 00340 eps = 2.22D-16 ! Set eps = zero here if you want it computed. 00341 00342 eps = 0.0D0 !!!!!!!!!!!!!!! 00343 gmin = 0.0D0 00344 gmax = 0.0D0 00345 00346 IF (eps <= zero) THEN 00347 eps = two**(-12) 00348 DO 00349 eps = eps / two 00350 x(1) = eps 00351 y(1) = one 00352 CALL daxpy ( 1, one, x, 1, y, 1 ) 00353 IF (y(1) <= one) EXIT 00354 END DO 00355 eps = eps * two 00356 END IF 00357 ! WRITE(*,*) 'computed epsilon is ',eps !!!!!!!!!!!!!!! 00358 00359 ! ------------------------------------------------------------------ 00360 ! Print heading and initialize. 00361 ! ------------------------------------------------------------------ 00362 IF (nout > 0) THEN 00363 WRITE(nout, 1000) enter, n, checka, precon, itnlim, rtol, shift 00364 END IF 00365 istop = 0 00366 itn = 0 00367 anorm = zero 00368 acond = zero 00369 rnorm = zero 00370 ynorm = zero 00371 CALL dload2( n, zero, x ) 00372 00373 ! ------------------------------------------------------------------ 00374 ! Set up y and v for the first Lanczos vector v1. 00375 ! y = beta1 P' v1, where P = C**(-1). 00376 ! v is really P' v1. 00377 ! ------------------------------------------------------------------ 00378 CALL dcopy ( n, b, 1, y , 1 ) ! y = b 00379 CALL dcopy ( n, b, 1, r1, 1 ) ! r1 = b 00380 IF ( precon ) CALL msolve( n, b, y ) 00381 beta1 = ddot ( n, b, 1, y, 1 ) 00382 00383 IF (beta1 < zero) THEN ! m must be indefinite. 00384 istop = 8 00385 GO TO 900 00386 END IF 00387 00388 IF (beta1 == zero) THEN ! b = 0 exactly. Stop with x = 0. 00389 istop = 0 00390 GO TO 900 00391 END IF 00392 00393 beta1 = SQRT( beta1 ) ! Normalize y to get v1 later. 00394 00395 ! ------------------------------------------------------------------ 00396 ! See if Msolve is symmetric. 00397 ! ------------------------------------------------------------------ 00398 IF (checka .AND. precon) THEN 00399 CALL msolve( n, y, r2 ) 00400 s = ddot ( n, y, 1, y, 1 ) 00401 t = ddot ( n,r1, 1,r2, 1 ) 00402 z = ABS( s - t ) 00403 epsa = (s + eps) * eps**0.33333D+0 00404 IF (z > epsa) THEN 00405 istop = 7 00406 GO TO 900 00407 END IF 00408 END IF 00409 00410 ! ------------------------------------------------------------------ 00411 ! See if Aprod is symmetric. 00412 ! ------------------------------------------------------------------ 00413 IF (checka) THEN 00414 CALL aprod ( n, y, w ) 00415 CALL aprod ( n, w, r2 ) 00416 s = ddot ( n, w, 1, w, 1 ) 00417 t = ddot ( n, y, 1,r2, 1 ) 00418 z = ABS( s - t ) 00419 epsa = (s + eps) * eps**0.33333D+0 00420 IF (z > epsa) THEN 00421 istop = 6 00422 GO TO 900 00423 END IF 00424 END IF 00425 00426 ! ------------------------------------------------------------------ 00427 ! Initialize other quantities. 00428 ! ------------------------------------------------------------------ 00429 oldb = zero 00430 beta = beta1 00431 dbar = zero 00432 epsln = zero 00433 qrnorm = beta1 00434 phibar = beta1 00435 rhs1 = beta1 00436 rhs2 = zero 00437 tnorm2 = zero 00438 ynorm2 = zero 00439 cs = - one 00440 sn = zero 00441 CALL dload2( n, zero, w ) ! w = 0 00442 CALL dload2( n, zero, w2 ) ! w2 = 0 00443 CALL dcopy ( n, r1, 1, r2, 1 ) ! r2 = r1 00444 00445 IF (debug) THEN 00446 WRITE(*,*) ' ' 00447 WRITE(*,*) 'b ', b 00448 WRITE(*,*) 'beta ', beta 00449 WRITE(*,*) ' ' 00450 END IF 00451 00452 ! ------------------------------------------------------------------ 00453 ! Main iteration loop. 00454 ! ------------------------------------------------------------------ 00455 DO 00456 itn = itn + 1 ! k = itn = 1 first time through 00457 IF (istop /= 0) EXIT 00458 00459 !----------------------------------------------------------------- 00460 ! Obtain quantities for the next Lanczos vector vk+1, k = 1, 2,... 00461 ! The general iteration is similar to the case k = 1 with v0 = 0: 00462 ! 00463 ! p1 = Operator * v1 - beta1 * v0, 00464 ! alpha1 = v1'p1, 00465 ! q2 = p2 - alpha1 * v1, 00466 ! beta2^2 = q2'q2, 00467 ! v2 = (1/beta2) q2. 00468 ! 00469 ! Again, y = betak P vk, where P = C**(-1). 00470 ! .... more description needed. 00471 !----------------------------------------------------------------- 00472 s = one / beta ! Normalize previous vector (in y). 00473 CALL dscal2( n, s, y, v ) ! v = vk if P = I 00474 00475 CALL aprod ( n, v, y ) 00476 CALL daxpy ( n, (- shift), v, 1, y, 1 ) 00477 IF (itn >= 2) THEN 00478 CALL daxpy ( n, (- beta/oldb), r1, 1, y, 1 ) 00479 END IF 00480 00481 alfa = ddot ( n, v, 1, y, 1 ) ! alphak 00482 00483 CALL daxpy ( n, (- alfa/beta), r2, 1, y, 1 ) 00484 CALL dcopy ( n, r2, 1, r1, 1 ) 00485 CALL dcopy ( n, y, 1, r2, 1 ) 00486 IF ( precon ) CALL msolve( n, r2, y ) 00487 00488 oldb = beta ! oldb = betak 00489 beta = ddot ( n, r2, 1, y, 1 ) ! beta = betak+1^2 00490 IF (beta < zero) THEN 00491 istop = 6 00492 WRITE(*,*) 'unfortunately beta < 0 ',beta,zero 00493 EXIT 00494 END IF 00495 00496 beta = SQRT( beta ) ! beta = betak+1 00497 tnorm2 = tnorm2 + alfa**2 + oldb**2 + beta**2 00498 00499 IF (itn == 1) THEN ! Initialize a few things. 00500 IF (beta/beta1 <= ten*eps) THEN ! beta2 = 0 or ~ 0. 00501 istop = -1 ! Terminate later. 00502 END IF 00503 !tnorm2 = alfa**2 00504 gmax = ABS( alfa ) ! alpha1 00505 gmin = gmax ! alpha1 00506 END IF 00507 00508 ! Apply previous rotation Qk-1 to get 00509 ! [deltak epslnk+1] = [cs sn][dbark 0 ] 00510 ! [gbar k dbar k+1] [sn -cs][alfak betak+1]. 00511 00512 oldeps = epsln 00513 delta = cs * dbar + sn * alfa ! delta1 = 0 deltak 00514 gbar = sn * dbar - cs * alfa ! gbar 1 = alfa1 gbar k 00515 epsln = sn * beta ! epsln2 = 0 epslnk+1 00516 dbar = - cs * beta ! dbar 2 = beta2 dbar k+1 00517 00518 ! Compute the next plane rotation Qk 00519 00520 agamma = SQRT( gbar**2 + beta**2 ) ! gammak 00521 cs = gbar / agamma ! ck 00522 sn = beta / agamma ! sk 00523 phi = cs * phibar ! phik 00524 phibar = sn * phibar ! phibark+1 00525 00526 IF (debug) THEN 00527 WRITE(*,*) ' ' 00528 WRITE(*,*) 'v ', v 00529 WRITE(*,*) 'alfa ', alfa 00530 WRITE(*,*) 'beta ', beta 00531 WRITE(*,*) 'gamma', agamma 00532 WRITE(*,*) 'delta', delta 00533 WRITE(*,*) 'gbar ', gbar 00534 WRITE(*,*) 'epsln', epsln 00535 WRITE(*,*) 'dbar ', dbar 00536 WRITE(*,*) 'phi ', phi 00537 WRITE(*,*) 'phiba', phibar 00538 WRITE(*,*) ' ' 00539 END IF 00540 00541 ! Update x. 00542 00543 denom = one/agamma 00544 00545 DO i = 1, n 00546 w1(i) = w2(i) 00547 w2(i) = w(i) 00548 w(i) = ( v(i) - oldeps*w1(i) - delta*w2(i) ) * denom 00549 x(i) = x(i) + phi * w(i) 00550 END DO 00551 00552 ! Go round again. 00553 00554 gmax = MAX( gmax, agamma) 00555 gmin = MIN( gmin, agamma) 00556 z = rhs1 / agamma 00557 ynorm2 = z**2 + ynorm2 00558 rhs1 = rhs2 - delta * z 00559 rhs2 = - epsln * z 00560 00561 ! Estimate various norms and test for convergence. 00562 00563 anorm = SQRT( tnorm2 ) 00564 ynorm = SQRT( ynorm2 ) 00565 epsa = anorm * eps 00566 epsx = anorm * ynorm * eps 00567 epsr = anorm * ynorm * rtol 00568 diag = gbar 00569 IF (diag == zero) diag = epsa 00570 00571 qrnorm = phibar 00572 rnorm = qrnorm 00573 00574 ! Estimate cond(A). 00575 ! In this version we look at the diagonals of R in the 00576 ! factorization of the lower Hessenberg matrix, Q * H = R, 00577 ! where H is the tridiagonal matrix from Lanczos with one 00578 ! extra row, beta(k+1) e_k^T. 00579 00580 acond = gmax / gmin 00581 00582 ! See if any of the stopping criteria are satisfied. 00583 ! In rare cases, istop is already -1 from above (Abar = const*I). 00584 00585 IF (istop == 0) THEN 00586 IF (itn >= itnlim ) istop = 5 00587 IF (acond >= 0.1D+0/eps) istop = 4 00588 IF (epsx >= beta1 ) istop = 3 00589 IF (qrnorm <= epsx ) istop = 2 00590 IF (qrnorm <= epsr ) istop = 1 00591 END IF 00592 00593 00594 ! See if it is time to print something. 00595 00596 IF (nout > 0) THEN 00597 prnt = .false. 00598 IF (n <= 40 ) prnt = .true. 00599 IF (itn <= 10 ) prnt = .true. 00600 IF (itn >= itnlim - 10) prnt = .true. 00601 IF (MOD(itn,10) == 0) prnt = .true. 00602 IF (qrnorm <= ten * epsx) prnt = .true. 00603 IF (qrnorm <= ten * epsr) prnt = .true. 00604 IF (acond >= 1.0D-2/eps ) prnt = .true. 00605 IF (istop /= 0 ) prnt = .true. 00606 00607 IF ( prnt ) THEN 00608 IF ( itn == 1) WRITE(nout, 1200) 00609 WRITE(nout, 1300) itn, x(1), qrnorm, anorm, acond 00610 IF (MOD(itn,10) == 0) WRITE(nout, 1500) 00611 END IF 00612 END IF 00613 00614 END DO 00615 ! ------------------------------------------------------------------ 00616 ! End of main iteration loop. 00617 ! ------------------------------------------------------------------ 00618 00619 ! Display final status. 00620 00621 900 IF (nout > 0) THEN 00622 WRITE(nout, 2000) EXIT, istop, itn, EXIT, anorm, acond, & 00623 EXIT, rnorm, ynorm 00624 WRITE(nout, 3000) EXIT, msg(istop) 00625 END IF 00626 00627 RETURN 00628 00629 00630 1000 FORMAT(// 1P, a, 5X, 'Solution of symmetric Ax = b' & 00631 / ' n =', i7, 5X, 'checkA =', l4, 12X, 'precon =', l4 & 00632 / ' itnlim =', i7, 5X, 'rtol =', e11.2, 5X, 'shift =', e23.14) 00633 1200 FORMAT(// 5X, 'itn', 8X, 'x(1)', 10X, & 00634 'norm(r)', 3X, 'norm(A)', 3X, 'cond(A)') 00635 1300 FORMAT(1P, i8, e19.10, 3E10.2) 00636 1500 FORMAT(1X) 00637 2000 FORMAT(/ 1P, a, 5X, 'istop =', i3, 14X, 'itn =', i8 & 00638 / a, 5X, 'Anorm =', e12.4, 5X, 'Acond =', e12.4 & 00639 / a, 5X, 'rnorm =', e12.4, 5X, 'ynorm =', e12.4) 00640 3000 FORMAT( a, 5X, a ) 00641 00642 END SUBROUTINE minres
1.8.0
