*+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++ * File minrestest.f * * These routines are for testing MINRES. * * 15 Jul 2003: Derived from symmlqtest.f. *+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++ subroutine Aprod ( n, x, y ) implicit double precision (a-h,o-z) double precision x(n), y(n) * ------------------------------------------------------------------ * Aprod computes y = A*x for some matrix A. * This is a simple example for testing MINRES. * ------------------------------------------------------------------ do 10 i = 1, n d = i !!!! * 1.1 d = d / n y(i) = d * x(i) 10 continue end ! subroutine Aprod *+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++ subroutine Msolve( n, x, y ) implicit double precision (a-h,o-z) double precision x(n), y(n) * ------------------------------------------------------------------ * Msolve solves M*y = x for some symmetric positive-definite * matrix M. * This is a simple example for testing MINRES. * shiftm will be the same as shift in MINRES. * * If pertbn = 0, the preconditioner will be exact, so * MINRES should require either one or two iterations, * depending on whether (A - shift*I) is positive definite or not. * * If pertbn is nonzero, somewhat more iterations will be required. * ------------------------------------------------------------------ intrinsic abs, mod common /mshift/ shiftm, pertm do 10 i = 1, n d = i !!!! * 1.1 d = d / n d = abs( d - shiftm ) if (mod(i,10) .eq. 0) d = d + pertm y(i) = x(i) / d 10 continue end ! subroutine Msolve *+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++ subroutine test ( n, precon, shift, pertbn ) implicit double precision (a-h,o-z) logical precon * ------------------------------------------------------------------ * test solves sets up and solves a system (A - shift * I)x = b, * using Aprod to define A and Msolve to define a preconditioner. * ------------------------------------------------------------------ common /mshift/ shiftm, pertm intrinsic abs external aprod, msolve logical checkA double precision b(100), r1(100), r2(100), $ v(100), w (100), w1(100), w2(100), $ x(100), y (100), xtrue(100) parameter ( one = 1.0d+0, two = 2.0d+0 ) shiftm = shift pertm = abs( pertbn ) nout = 6 write(nout, 1000) shift, pertbn * Set the true solution and the rhs * so that (A - shift*I) * xtrue = b. do 100 i = 1, n xtrue(i) = n + 1 - i 100 continue call aprod ( n, xtrue, b ) call daxpy ( n, (- shift), xtrue, 1, b, 1 ) * Set other parameters and solve. checkA = .true. itnlim = n * 2 rtol = 1.0d-12 call minres( n, b, r1, r2, v, w, w1, w2, x, y, $ Aprod, Msolve, checkA, precon, shift, $ nout , itnlim, rtol, $ istop, itn, Anorm, Acond, rnorm, ynorm ) * Compute the final residual, r1 = b - (A - shift*I)*x. call Aprod ( n, x, y ) call daxpy ( n, (- shift), x, 1, y, 1 ) do 400 i = 1, n r1(i) = b(i) - y(i) 400 continue r1norm = dnrm2 ( n, r1, 1 ) write(nout, 2000) r1norm * Print the solution and some clue about whether it is OK. write(nout, 2100) (j, x(j), j = 1, n) do 500 j = 1, n w(j) = x(j) - xtrue(j) 500 continue enorm = dnrm2 ( n, w, 1 ) / dnrm2 ( n, xtrue, 1 ) etol = 1.0d-5 if (enorm .le. etol) then write(nout, 3000) enorm else write(nout, 3100) enorm end if 1000 format(//' ------------------------------------------------------' $ / ' Test of MINRES.' $ / ' ------------------------------------------------------' $ / ' shift =', f12.4, 6x, 'pertbn =', f12.4) 2000 format(/ ' Final residual =', 1p, e8.1) 2100 format(/ ' Solution x' / 1p, 4(i6, e14.6)) 3000 format(/ ' MINRES appears to be successful.', $ ' Relative error in x =', 1p, e8.1) 3100 format(/ ' MINRES appears to have failed. ', $ ' Relative error in x =', 1p, e8.1) end ! subroutine test * ------------------------------------------------------------------ * Main Program. * ------------------------------------------------------------------ program main implicit double precision (a-h,o-z) logical normal, precon parameter ( zero = 0.0d+0, one = 1.0d+0 ) normal = .false. precon = .true. shift = 0.25 pertbn = 0.1 * Test the unlikely tiny cases that often trip us up. n = 1 call test( n, normal, zero , zero ) call test( n, normal, shift, zero ) n = 2 call test( n, normal, zero , zero ) call test( n, normal, shift, zero ) * Test small positive-definite and indefinite systems * without preconditioners. MINRES should take n iterations. n = 50 call test( n, normal, zero , zero ) call test( n, normal, shift, zero ) ! if (n .le. 5) stop ! WHILE TESTING * Test small positive-definite and indefinite systems with * exact preconditioners. MINRES should take 1 and 2 iterations * respectively. n = 1 call test( n, precon, zero , zero ) call test( n, precon, shift, zero ) n = 2 call test( n, precon, zero , zero ) call test( n, precon, shift, zero ) n = 50 call test( n, precon, zero , zero ) call test( n, precon, shift, zero ) * pertbn makes the preconditioners incorrect in n/10 entries. * MINRES should take about n/10 iterations. call test( n, precon, zero , pertbn ) call test( n, precon, shift, pertbn ) end ! program main