In this section we give an example of the layout of a subroutine according to the rules given above.
subroutine profun ( amat, intdia, intpro, n, nmat )
c ======================================================================
c
c
c programmer Guus Segal
c version 1.1 date 26-11-93 (New type of error messages)
c version 1.0 date 10-11-91
c
c
c
c copyright (c) 1991,1993 "Ingenieursbureau SEPRA"
c permission to copy or distribute this software or documentation
c in hard copy or soft copy granted only by written licence
c obtained from "Ingenieursbureau SEPRA".
c all rights reserved. no part of this publication may be reproduced,
c stored in a retrieval system ( e.g., in memory, disk, or core)
c or be transmitted by any means, electronic, mechanical, photocopy,
c recording, or otherwise, without written permission from the
c publisher.
c
c
c ********************************************************************
c
c DESCRIPTION
c
c SEPRAN subroutine to compute the L U decomposition of a
c non-symmetric real matrix.
c The matrix is supposed to be stored as a profile matrix.
c
c
c The result of the decomposition may depend on the accuracy of the
c computer and the condition of the matrix.
c There is no pivoting applied in this subroutine.
c
c For a 32-bits machine double precision is almost always necessary in
c order to guarantee accuracy.
c ********************************************************************
c
c KEYWORDS
c
c linear_solver
c lu_decomposition
c profile_storage
c ********************************************************************
c
c INPUT / OUTPUT PARAMETERS
c
implicit none
integer n, intdia(n), intpro(n), nmat
double precision amat(nmat)
c amat i/o In this real array the parts of the matrix to be
c decomposed are stored at input. At output the decomposition
c overwrites the input
c For the matrix a profile storage is assumed as described
c in the SEPRAN Programmer's Guide 13.6
c intdia i Integer array of length n containing the positions of the
c diagonal elements in amat, i.e. A = amat ( intdia(i) )
c ii
c intpro i Integer array of length n containing the column numbers
c of the at most left elements in a row in amat
c n i Number of rows of the matrix (n x n )
c nmat i Size of the matrix
c ********************************************************************
c
c COMMON BLOCKS
double precision epsmac, sqreps, rinfin
common /cmcdpr/ epsmac, sqreps, rinfin
c
c /cmcdpr/
c Machine dependent constants.
c
c epsmac i Machine accuracy
c sqreps i Square root of machine accuracy
c rinfin i Approximation of infinity
c ********************************************************************
c
c LOCAL PARAMETERS
c
integer icolpi, icol, ii, i, j, jj, ij, ji
double precision rowsum
c i Counting variable
c icol Max (pi, pj)
c icolpi Column number of at most left non-zero element in row i
c ii Position of i th diagonal element in array amat
c ij Position of element ij in array amat
c j Counting variable
c ji Position of element ji in array amat
c jj Position of j th diagonal element in array amat
c rowsum Row norm of matrix
c *******************************************************************
c
c I/O
c
c none
c ********************************************************************
c
c SUBROUTINES CALLED
double precision dasum, ddot
c DASUM (SASUM) Sum of magnitudes of vector components
c DDOT (CDOT) Inner product of two vectors
c ERCLOS Resets old name of previous subroutine of higher level
c EROPEN Produces concatenated name of local subroutine
c ERREAL Put real in error message
c ERRINT Put integer in error message
c ERRSUB Error messages
c INSTOP Stop the program
c ********************************************************************
c
c ERROR MESSAGES
c
c 104 matrix is singular
c 1820 The matrix found in the solver is completely zero
c ********************************************************************
c
c PSEUDO CODE
c
c
c METHOD
c
c The L U decomposition of the matrix A is based on the following
c formulas:
c
c n
c a = sum l u l = 1
c ij k=1 ik kj kk
c
c this implies:
c
c j-1
c u = a - sum l u ( j=P (1) i-1 )
c ji ji k=max(P ,P ) jk ki i
c i j
c
c
c j-1
c l = ( a - sum l u ) / d ( j=P (1) i-1 )
c ij ij k=max(P ,P ) ik kj jj i
c i j
c
c
c i-1
c u = d = a - sum l u
c ii ii ii k=P ik ki
c i
c
c Here P denotes the column number of the first non-zero element in row i
c i
c
c
c The decomposition is stopped if |d | < ||A|| * epsmac
c ii
c
c For the matrix norm we use the maximum row norm, i.e.
c
c ||A|| = max sum |a |
c i j ij
c
c
c Pseudo code
c
c Compute row norm of matrix and multiply by machine accuracy
c For i = 1 (1) n
c For j = P (1) i-1
c i
c
c P = max ( P , P )
c i j
c
c If ( P < j ) then
c
c adapt u with inner product
c ji
c
c adapt l with inner product
c ij
c
c l = l / d
c ij ij ii
c
c adapt d with inner product
c ii
c
c if (d < ||A|| epsmac ) then
c ii
c
c Give error message and stop the program
c ======================================================================
c
c --- compute row norm of part of the matrix in order to be able to check
c the stopping criterion for the pivot process
c
rowsum = 0d0
do 100 i = 1, n
rowsum = max ( rowsum, dasum (2*i-intpro(i)+1,
+ amat(intdia(i)-i-intpro(i)), 1) )
100 continue
if ( rowsum.eq.0d0 ) then
call eropen ( 'profun' )
call errsub ( 1820, 0, 0, 0 )
call erclos ( 'profun' )
call instop
end if
rowsum = rowsum*epsmac
c --- decomposition algorithm
c Row loop
do 150 i = 1, n
ii = intdia(i)
icolpi = intpro(i)
if ( icolpi.lt.i ) then
c --- pi < i
c Loop over columns
do 110 j = icolpi, i-1
icol = max ( icolpi, intpro(j) )
jj = intdia(j)
ij = ii-i+j
if ( icol.lt.j ) then
c --- p < j
ji=ii+i-j
amat(ji) = amat(ji) - ddot ( j-icol, amat(jj-j+icol),
+ 1, amat(ji+1), -1 )
amat(ij) = amat(ij) - ddot ( j-icol, amat(ii-i+icol),
+ 1, amat(jj+1), -1 )
end if
amat(ij) = amat(ij)/amat(jj)
110 continue
c --- Diagonal element
amat(ii) = amat(ii) - ddot ( i-icolpi, amat(ii-i+icolpi),
+ 1, amat(ii+1), -1 )
end if
c --- Check pivot element:
if ( abs(amat(ii)).lt.rowsum ) then
call eropen ( 'profun' )
call errint ( i, 1 )
call erreal ( amat(ii), 1 )
call erreal ( rowsum, 2 )
call errsub ( 104, 1, 2, 0 )
call erclos ( 'profun' )
call instop
end if
150 continue
c --- end of decomposition
end