subroutine dbiplr(a,b,m,bda,bdb,n,f,idf,
+ alpha,beta,iflag,w,lw)
c
integer m,n,idf,iflag,lw
double precision a,b,alpha,beta
double precision bda(*),bdb(*),f(idf,*),w(*)
c
c this subroutine solves the generalized biharmonic
c equation
c
c
c 1 1 2 1 1 2
c - ( d r d + - d ) - ( d r d + - d ) u(r,t) +
c r r r r t r r r r t
c
c 1 1 2
c alpha ( - ( d r d + - d ) u(r,t) + beta u(r,t) = f(r,t)
c r r r r t
c
c
c in polar coordinates. the boundary conditions are of the first
c kind, u and the first derivative of u with respect to r must
c be specified on the boundary.
c r is radial coordinate, while t denotes the theta coordinate.
c the equation is solved in a region between two concentric
c circles. the radius of the inner circle being a and that
c of the outer circle being b.
c the special case when the region
c is a disk is accepted when a is set to zero.
c
c the solution time is proportional to m*n*log(n).
c
c this code is written by petter bjorstad, the author
c acknowledges paul n swarztrauber for the use of his
c fast fourier transform package and the linpack project
c for the use of linear equation solvers.
c
c this is version 2.0 , double precision. january 1984.
c
c any comments, questions or suggestions are most welcome
c and should be directed to :
c
c petter bjorstad
c
c veritas research
c p.o. box 300
c n-1322 hovik
c norway
c
c description of input variables.
c
c a,b defines the domain where the equation is defined.
c see explanation above.
c note the special meaning if a is set to zero.
c
c m the number of interior gridpoints in the
c r direction. the interval a .le. x .le. b is divide
c in m+1 panels each of width (b-a)/(m+1).
c m must be at least 3.
c
c n the number of gridpoints in the theta direction.
c the angle 2 pi (6.28...) is divided into n panels
c each of width 2*pi/n. n should be at least 3.
c the method is more efficient if n is a product of
c small primes.
c
c alpha the constant alpha in the equation.
c
c beta the constant beta in the equation.
c
c bda array of dimension at least n.
c the dimension must be at least m if a is zero.
c if a is non-zero, then
c bda(j) contains the derivative with respect
c to r of the function at the inner boundary.
c (at r = a and t = j*2*pi/n ,j=1,2,..n.)
c if the parameter a is zero then this array
c of length at least m is used as an extra
c workspace.
c
c bdb array of dimension at least n.
c bdb(j) contains the derivative with respect
c to r of the function at the outer boundary.
c (at r = b and t = j*2*pi/n ,j=1,2,..n.)
c
c f array of dimension at least ( m+2,n ) .
c f(i,j) must contain the values of the right
c hand side function f(r,t) evaluated at the points
c (r ,t ) , where
c i j
c r = a+(i-1)*(b-a)/(m+1) , i=2,3,.. m+1.
c i
c t = j*2*pi/n , j=1,2,.. n.
c j
c f(1,j) contains the boundary values at the
c inner boundary. (r =a , t = t )
c j
c if a is zero then f(1,j) should contain the value
c of the function f(0,t). (the element f(1,1) will
c be read by the code).
c
c f(m+2,j) contains the boundary values at the
c outer boundary. (r = b , t = t )
c j
c
c idf rowdimension of f as declared in the
c calling program.
c
c iflag = 1 if the discrete system is positive definite.
c this is always the case if alpha is less or equal
c zero and beta is greater or equal zero. this option
c is somewhat faster than using iflag = 2.
c = 2 this option should handle all nonsingular discrete
c versions of the equations.
c
c description of output variables.
c
c f contains the solution u of the discrete
c system in all gridpoints (r ,t ) , i=1,2,.. m+2.
c i j , j=1,2,.. n.
c
c iflag error indicator.
c = 0 normal return.
c =-1 n and/or m was less than 3.
c =-2 a is greater or equal b or a is negative.
c =-3 idf is less than m+2 or lw is too small.
c =-4 linpack failure in cholesky factorization.
c this can happen if the routine is called with
c iflag=1 and the linear system is indefinite.
c try calling with iflag=2.
c this error should never occur if alpha is less or
c equal zero and beta is greater or equal zero.
c =-5 linpack detected a computationally singular
c system using gaussian elimination with pivoting.
c
c description of workspace.
c
c lw the length of the user supplied workspace w.
c lw depends on iflag in the following way:
c
c iflag lw must be at least
c
c 1 2*m+n+max(8*m+4,2*n+15)
c 2 2*m+n+max(13*m,2*n+15)
c
c w a one dimensional array containing at least lw
c elements. w can be used for other purposes
c between calls to dbiplr.
c
c subroutines and functions
c needed by this program.
c
c 1. biharmonic: dbipl
c 2. fourier
c transform : drffti, drfti1, drfftf, drftf1,
c dradf2, dradf3, dradf4, dradf5,
c dradfg, drfftb, drftb1, dradb2,
c dradb3, dradb4, dradb5, dradbg.
c
c 3. linpack : dpbfa, dpbsl, dgbfa, dgbsl.
c
c 4. blas : dscal, dcopy, ddot, daxpy, idamax
c
c 5. fortran : min, max, mod, cos, atan, sqrt
c sin
c
c
c local.
c
c biharmonic: dbipl
c fortran : max
c
integer i1,i2,i3,i4,i5,i6,i7,i8
c
c check for input errors.
c
i1 = 2*m + n + max(8*m+4, 2*n+15)
i2 = 2*m + n + max(13*m , 2*n+15)
c
if(n .lt. 3 .or. m .lt. 3) iflag = -1
if(a .ge. b) iflag = -2
if(a .lt. 0.0d0) iflag = -2
if(idf .lt. m+2) iflag = -3
if(iflag .eq. 1 .and. lw .lt. i1) iflag = -3
if(iflag .eq. 2 .and. lw .lt. i2) iflag = -3
if(iflag .lt. 0) go to 100
c
c assign workspace.
c
i1 = 1
i2 = m+i1
i3 = m+i2
i4 = n+i3
i5 = m+i4
i6 = m+i5
i7 = m+i6
i8 = m+i7
c
call dbipl(a,b,m,bda,bdb,n,f,idf,alpha,beta,iflag,
+ w(i1),w(i2),w(i3),w(i4),w(i5),w(i6),
+ w(i7),w(i8))
c
if(iflag.ge.0) return
100 write(6,1) iflag
1 format(1x,'error return from dbiplr , iflag= ',i4)
end