*+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++ * File symmlqtest.f * * These routines are for testing SYMMLQ. * * 04 Sep 1991: logical checkA added to parameter list. *+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++ 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 SYMMLQ. * ------------------------------------------------------------------ do 10 i = 1, n d = i * 1.1 d = d / n y(i) = d * x(i) 10 continue * end of Aprod end 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 SYMMLQ. * shiftm will be the same as shift in SYMMLQ. * * If pertbn = 0, the preconditioner will be exact, so * SYMMLQ 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 of Msolve end subroutine test ( n, goodb, precon, shift, pertbn ) implicit double precision (a-h,o-z) logical goodb, 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) double precision 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.0E-12 call symmlq( n, b, r1, r2, v, w, x, y, $ aprod, msolve, checkA, goodb, 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.0e-5 if (enorm .le. etol) then write(nout, 3000) enorm else write(nout, 3100) enorm end if 1000 format(//' ------------------------------------------------------' $ / ' Test of SYMMLQ.' $ / ' ------------------------------------------------------' $ / ' 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(/ ' SYMMLQ appears to be successful.', $ ' Relative error in x =', 1p, e8.1) 3100 format(/ ' SYMMLQ appears to have failed. ', $ ' Relative error in x =', 1p, e8.1) * end of test end * ------------------------------------------------------------------ * Main Program. * ------------------------------------------------------------------ implicit double precision (a-h,o-z) logical goodb, normal, precon parameter ( zero = 0.0d+0, one = 1.0d+0 ) normal = .false. precon = .true. shift = 0.25 pertbn = 0.1 do 100 ntest = 1, 2 if (ntest .eq. 1) goodb = .false. if (ntest .eq. 2) goodb = .true. * Test the unlikely tiny cases that often trip us up. call test( 1, goodb, normal, zero , zero ) call test( 2, goodb, normal, zero , zero ) call test( 1, goodb, precon, zero , zero ) call test( 2, goodb, precon, zero , zero ) * Test small positive-definite and indefinite systems with * exact preconditioners. SYMMLQ should take 1 and 2 iterations * respectively. n = 5 call test( n, goodb, precon, zero , zero ) call test( n, goodb, precon, shift, zero ) * Test other combinations on larger systems. * With no preconditioning, SYMMLQ will take about n iterations. n = 50 call test( n, goodb, normal, zero , zero ) call test( n, goodb, normal, shift, zero ) * PERTBN makes the preconditioners incorrect in n/10 entries. * SYMMLQ should take about n/10 iterations. call test( n, goodb, precon, zero , pertbn ) call test( n, goodb, precon, shift, pertbn ) 100 continue * end of Main Program for testing SYMMLQ end