***** File LUMODTEST FORTRAN * * Main lutest setmat luchek random ddrand * * ++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++ program LUmodTest implicit double precision (a-h,o-z) parameter ( maxmod = 20, $ nrowb = 2*maxmod, ncolb = 2*maxmod ) double precision B(nrowb,ncolb), D(maxmod), x(maxmod) double precision C(maxmod,maxmod) double precision L(maxmod*maxmod), U(maxmod*(maxmod+1)/2) double precision v(maxmod), w(maxmod) double precision y(maxmod), z(maxmod) * ------------------------------------------------------------------ * Test program for LUMOD. * * 25 Apr 1990: First version, derived from test program for QRMOD. * LUTEST is asked to do mtest tests, each of which maintains the * factors of C as various modes of LUMOD are exercised. * C builds up to a maximum dimension of maxmod. * To make the test more demanding, increase maxmod and mtest. * * 22 Aug 1991: More statistics are now accumulated by LUTEST. * 29 Jun 1999: dzero replaced by dload. * 29 Jun 1999: Column replacement now done without reordering. * 30 Jun 1999: Row and column delete swaps with last row and col. * ------------------------------------------------------------------ integer seeds(3) double precision r1, r2, densty * Initialize mtest and the random number sequence. mtest = 2 r1 = -10.0d+0 r2 = 10.0d+0 densty = 0.2d+0 seeds(1) = 123456 seeds(2) = 234567 seeds(3) = 345678 * Test square systems. m = maxmod n = maxmod call lutest( m, n, nrowb, ncolb, mtest, $ r1, r2, densty, seeds, $ B, C, D, L, U, v, w, x, y, z ) * End of main program for LUTEST end * ++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++ subroutine lutest( m, n, nrowb, ncolb, mtest, $ r1, r2, densty, seeds, $ B, C, D, L, U, v, w, x, y, z ) implicit double precision (a-h,o-z) double precision B(nrowb,ncolb), C(m,m), D(m), x(m), y(m), z(m) double precision L(*), U(*) double precision v(m), w(m) integer seeds(3) double precision r1, r2, densty * ------------------------------------------------------------------ * LUTEST is a test program for LUMOD. * ------------------------------------------------------------------ double precision stats(7) double precision zero, one, big data zero /0.0d+0/, one /1.0d+0/, big /1.0d+30/ * Initialize nout, stats(*). * nout = file number for printing stats. * stats(1) = Cmax = max |Cij| * stats(2) = Lmax = max |Lij| * stats(3) = Umax = max |Uij| * stats(4) = Umin = min |Ujj| * stats(5) = Fmax = max |L*C - U| * stats(6) = smax = |Cx - b| * stats(7) = tmax = |C'y - b| nout = 9 stats(1) = zero stats(2) = zero stats(3) = zero stats(4) = big stats(5) = zero stats(6) = zero stats(7) = zero * Run the specified number of tests. do 600 ntest = 1, mtest write(nout, '(// a / a, i4, a, i6 / a)') $ ' ------------------------------------------', $ ' LUTEST', ntest, '. m =', m, $ ' ------------------------------------------' call setmat( m, n, nrowb, B(1, 1), D, v, w, $ r1, r2, densty, seeds, y ) call setmat( m, n, nrowb, B(1,n+1), D, v, w, $ r1, r2, densty, seeds, y ) call setmat( m, n, nrowb, B(m+1,1), D, v, w, $ r1, r2, densty, seeds, y ) * --------------------------------------------------------------- * Test mode 1, adding 1 row and a column at a time. * --------------------------------------------------------------- write(nout, '(// a)') ' LUmod (mode 1):' do 200 nn = 1, m do 180 k = 1, nn y(k) = B(nn,k) z(k) = B(k,nn) 180 continue call dcopy ( nn, y, 1, C(nn,1), m ) call dcopy ( nn, z, 1, C(1,nn), 1 ) call LUmod ( 1, m, nn, krow, kcol, L, U, y, z, w ) write(nout, '(/ a, i6, a)') ' Row and column', nn, $ ' added' call luchek( m, n, nn, C, L, U, $ nout, stats, v, w, x, y ) 200 continue * --------------------------------------------------------------- * Replace columns of the LU one by one. * Old ones are replaced by new ones with no reordering. * --------------------------------------------------------------- write(nout, '(// a)') ' LUmod (mode 2):' do 300 j = 1, n jrep = j jnew = n + j do 250 i = 1, m z(i) = B(i,jnew) 250 continue call dcopy ( m, z, 1, C(1,jrep), 1 ) call LUmod ( 2, m, m, krow, jrep, L, U, y, z, w ) write(nout, '(/ a, i6, a, i6)') ' Column', jrep, $ ' replaced by column', jnew call luchek( m, n, m, C, L, U, $ nout, stats, v, w, x, y ) 300 continue * --------------------------------------------------------------- * Replace rows of the LU one by one. * Old ones are replaced by new ones with no reordering. * --------------------------------------------------------------- write(nout, '(// a)') ' LUmod (mode 3):' do 400 i = 1, m irep = i inew = m + i do 350 j = 1, n y(j) = B(inew,j) 350 continue call dcopy ( n, y, 1, C(irep,1), m ) call LUmod ( 3, m, m, irep, kcol, L, U, y, z, w ) write(nout, '(/ a, i6, a, i6)') ' Row', irep, $ ' replaced by row', inew call luchek( m, n, m, C, L, U, $ nout, stats, v, w, x, y ) 400 continue * --------------------------------------------------------------- * Delete a specified row and column of the LU one by one, * replacing them by the last row and column. * --------------------------------------------------------------- write(nout, '(// a)') ' LUmod (mode 4):' nn = m do 500 kdel = 1, m - 1 krow = max( nn/2, 1 ) kcol = krow call dcopy ( nn , C(1,nn), 1, C(1,kcol), 1 ) call dcopy ( nn-1, C(nn,1), m, C(krow,1), m ) call LUmod ( 4, m, nn, krow, kcol, L, U, y, z, w ) nn = nn - 1 write(nout, '(/ a, i6, a, i6, a)') ' Row', krow, $ ' and column', kcol, ' deleted' call luchek( m, n, nn, C, L, U, $ nout, stats, v, w, x, y ) 500 continue 600 continue write(nout, '(/ 2a / a, i4, 2a /)') $ ' -----------------------------------------', $ '------------------------------------------', $ ' LUTEST statistics:' write(nout, 1000) write(nout, 1100) mtest, stats write(nout, '(2a)') $ ' -----------------------------------------', $ '------------------------------------------' emax = max( stats(5), stats(6), stats(7) ) if (emax .le. 1.0d-6) then write(nout, '(/ a, 1p, e11.1)') $ ' Everything seems OK. Max error =', emax else write(nout, '(/ a, 1p, e11.1)') $ ' There may be a bug. Max error =', emax end if return 1000 format(' Tests max|Cij| max|Lij| max|Uij| min|Ujj|', $ ' max|LC-U| max|Cx-b| max|C''y-b|') 1100 format(i6, 1p, 7e11.1) * end of lutest end * ++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++ subroutine setmat( m, n, nrowb, B, D, v, w, $ r1, r2, densty, seeds, y ) implicit double precision (a-h,o-z) integer n, seeds(3) double precision B(nrowb,n), D(n), v(m), w(n), y(n) double precision r1, r2, densty * ------------------------------------------------------------------ * setmat sets B = D + vw', * where D = diag(d) and d, v, w are random vectors. * ------------------------------------------------------------------ dr1 = (r2 - r1) + 1.0d+0 dr2 = dr1 + 1.0d+0 ddnsty = 1.0d+0 call random( n, D, dr1, dr2, ddnsty, seeds, y ) call random( n, v, r1, r2, densty, seeds, y ) call random( n, w, r1, r2, densty, seeds, y ) do 50 j = 1, n wj = w(j) do 40 i = 1, m B(i,j) = v(i) * wj 40 continue B(j,j) = B(j,j) + D(j) 50 continue * end of setmat end * ++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++ subroutine luchek( m, n, nn, C, L, U, $ nout, stats, v, w, x, y ) implicit double precision (a-h,o-z) double precision C(m,n), x(n), y(m) double precision L(*), U(*) double precision v(m), w(n), stats(7) * ------------------------------------------------------------------ * LUCHEK computes L*C and compares it with U. * ------------------------------------------------------------------ character*1 mark double precision Cmax, Lmax, Umax, Umin, Fmax, smax, tmax double precision zero, one, big data zero /0.0d+0/, one /1.0d+0/, big /1.0d+30/ * ------------------------------------------------------------------ * Get Cmax, Lmax, Umax, Umin. * Prevent divide by zero in Usolve -- make sure U(k,k) .ne. 0 * ------------------------------------------------------------------ Cmax = zero Lmax = zero Umax = zero Umin = big ll = 1 lu = 1 incl = m incu = m numu = nn do 10 k = 1, nn icmax = idamax( nn , C(1,k), 1 ) ilmax = idamax( nn , L(ll), 1 ) iumax = idamax( numu, U(lu), 1 ) Cmax = max( Cmax, abs(C(icmax,k)) ) Lmax = max( Lmax, abs(L(ilmax)) ) Umax = max( Umax, abs(U(iumax)) ) Umin = min( Umin, abs(U(lu)) ) if (U(lu) .eq. zero) U(lu) = 1.0d-16 ll = ll + incl lu = lu + incu incu = incu - 1 numu = numu - 1 10 continue * ------------------------------------------------------------------ * Check L*C = U by columns. * ------------------------------------------------------------------ Fmax = zero do 100 k = 1, nn * Set v = L * k-th column of C. call Lprod ( 1, m, nn, L, C(1,k), v ) * Find the maximum deviation from the corresponding column of U. call dload ( nn, zero, w, 1 ) lu = k incu = m - 1 do 50 i = 1, k w(i) = U(lu) lu = lu + incu incu = incu - 1 50 continue do 60 i = 1, nn err = abs( w(i) - v(i) ) Fmax = max( Fmax, err ) 60 continue 100 continue * ------------------------------------------------------------------ * Check solution of C*x = v. * ------------------------------------------------------------------ x(1) = one do 300 j = 2, nn x(j) = - x(j-1) / 2.0d+0 300 continue call dload ( nn, zero, v, 1 ) do 320 k = 1, nn call daxpy ( nn, x(k), C(1,k), 1, v, 1 ) 320 continue * Save v in y. * Solve C*x = v = y, * then compute the residual y - C*x. call dcopy ( nn, v, 1, y, 1 ) call Lprod ( 1, m, nn, L, v, x ) call Usolve( 1, m, nn, U, x ) do 350 k = 1, nn call daxpy ( nn, -x(k), C(1,k), 1, y, 1 ) 350 continue imax = idamax( nn, y, 1 ) smax = abs ( y(imax) ) * ------------------------------------------------------------------ * Check solution of C(t)*y = w. * ------------------------------------------------------------------ y(1) = one do 400 i = 2, nn y(i) = - y(i-1) / 2.0d+0 400 continue do 420 k = 1, nn w(k) = ddot ( nn, C(1,k), 1, y, 1 ) 420 continue * Save w in x. * Solve C(t)*y = w = x, * then compute the residual x - C(t)*y. call dcopy ( nn, w, 1, x, 1 ) call Usolve( 2, m, nn, U, w ) call Lprod ( 2, m, nn, L, w, y ) do 450 k = 1, nn x(k) = x(k) - ddot ( nn, C(1,k), 1, y, 1 ) 450 continue imax = idamax( nn, x, 1 ) tmax = abs ( x(imax) ) mark = ' ' if ( max(smax,tmax) .gt. 1.0d-4 ) mark = '*' * Print stats for this matrix * and save max or min so far. write(nout, 1000) write(nout, 1100) $ nn, Cmax, Lmax, Umax, Umin, Fmax, smax, tmax, mark stats(1) = max( stats(1), Cmax ) stats(2) = max( stats(2), Lmax ) stats(3) = max( stats(3), Umax ) stats(4) = min( stats(4), Umin ) stats(5) = max( stats(5), Fmax ) stats(6) = max( stats(6), smax ) stats(7) = max( stats(7), tmax ) return 1000 format(' n max|Cij| max|Lij| max|Uij| min|Ujj|', $ ' max|LC-U| max|Cx-b| max|C''y-b|') 1100 format(i6, 1p, 7e11.1, 1x, a) * end of luchek end * ++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++ subroutine random( n, x, r1, r2, densty, seeds, y ) implicit none integer n, seeds(3) double precision x(n), y(n) double precision r1, r2, densty * ------------------------------------------------------------------ * random generates a vector x(*) of random numbers * in the range (r1, r2) with specified density. * seeds(*) must be initialized before the first call. * ------------------------------------------------------------------ integer i call ddrand( n, x, 1, seeds ) call ddrand( n, y, 1, seeds ) do i = 1, n if (y(i) .lt. densty) then x(i) = r1 + (r2 - r1) * x(i) else x(i) = 0.0d+0 end if end do * end of random end * ++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++ subroutine ddrand( n, x, incx, seeds ) implicit none integer n, incx, seeds(3) double precision x(*) * ------------------------------------------------------------------ * ddrand fills a vector x with uniformly distributed random numbers * in the interval (0, 1) using a method due to Wichman and Hill. * * seeds(1:3) should be set to integer values * between 1 and 30000 before the first entry. * * Integer arithmetic up to 30323 is required. * * Blatantly copied from Wichman and Hill 19-January-1987. * 14-Feb-94. Original version. * 30 Jun 1999. seeds stored in an array. * 30 Jun 1999. This version of ddrand. * ------------------------------------------------------------------ integer ix if (n .lt. 1) return do ix = 1, 1+(n-1)*incx, incx seeds(1) = 171*mod(seeds(1), 177) - 2*(seeds(1)/177) seeds(2) = 172*mod(seeds(2), 176) - 35*(seeds(2)/176) seeds(3) = 170*mod(seeds(3), 178) - 63*(seeds(3)/178) if (seeds(1) .lt. 0) seeds(1) = seeds(1) + 30269 if (seeds(2) .lt. 0) seeds(2) = seeds(2) + 30307 if (seeds(3) .lt. 0) seeds(3) = seeds(3) + 30323 x(ix) = mod( real(seeds(1))/30269.0 + $ real(seeds(2))/30307.0 + $ real(seeds(3))/30323.0, 1.0 ) end do * end of ddrand end