subroutine leg(hoyre,coeff,xjac,alf,bet,flag,N,x1,x2,tid1,tid2,tid3,tid4) ! ! ! Solves the generalized biharmonic equation using Shen's ! new Legendre basis. The solver used in this routine is ! symmetric, while Shen's propsed solver is non-symmetric. ! ! use Interface, ONLY:leg_hele,basistrans,leg_disk,leg_eig,leg_sepfakt, & leg_pfirst, leg_capmatr,dpotrf, sptoph2,phtosp2,legcoef,dpbtrf ! use Bi_Parameter ! implicit none ! integer, intent(in) :: flag,N double precision, intent(in) :: alf,bet double precision, intent(in) :: xjac(0:2*N+3) double precision, intent(in) :: x1,x2 double precision, intent(inout) :: hoyre(0:2*N+3,0:2*N+3) double precision, intent(out) :: coeff(0:2*N+3,0:2*N+3) ! ! integer i,ii,l,info,odde ! integer, save :: oldn=0 ! double precision k ! ! for timing purposes only ! integer count,ocount,count_rate,count_max double precision, intent(out) :: tid1,tid2,tid3,tid4 ! ! Variables that are saved for subsequent right hand sides. ! double precision, dimension(:,:),allocatable, save :: Id1 double precision, dimension(:,:),allocatable, save :: Id2 double precision, dimension(:,:),allocatable, save :: Id1trans double precision, dimension(:,:),allocatable, save :: Id2trans double precision, dimension(:,:),allocatable, save :: Ja1 double precision, dimension(:,:),allocatable, save :: Ja2 double precision, dimension(:,:),allocatable, save :: Fa1 double precision, dimension(:,:),allocatable, save :: Fa2 double precision, dimension(:,:),allocatable, save :: Sum1 double precision, dimension(:,:),allocatable, save :: Sum2 double precision, dimension(:,:,:),allocatable, save :: Stor1Sep double precision, dimension(:,:,:),allocatable, save :: Stor2Sep double precision, dimension(:,:,:),allocatable, save :: Stor3Sep double precision, dimension(:,:,:),allocatable, save :: Stor4Sep double precision, dimension(:,:),allocatable, save :: Cap1 double precision, dimension(:,:),allocatable, save :: Cap2 double precision, dimension(:,:),allocatable, save :: Cap3 double precision, dimension(:,:),allocatable, save :: Cap4 double precision, dimension(:),allocatable, save :: De1 double precision, dimension(:),allocatable, save :: De2 double precision, dimension(:),allocatable, save :: De3 double precision, dimension(:),allocatable, save :: De4 double precision, dimension(:),allocatable, save :: d double precision, dimension(:),allocatable, save :: r1 double precision, dimension(:),allocatable, save :: r2 double precision, dimension(:),allocatable, save :: r3 double precision, dimension(:),allocatable, save :: s2 double precision, dimension(:),allocatable, save :: s3 double precision, dimension(:),allocatable, save :: cap1_vekt double precision, dimension(:),allocatable, save :: cap2_vekt double precision, dimension(:),allocatable, save :: cap3_vekt double precision, dimension(:),allocatable, save :: cap4_vekt integer , dimension(:),allocatable, save :: info1 integer , dimension(:),allocatable, save :: info2 integer , dimension(:),allocatable, save :: info3 integer , dimension(:),allocatable, save :: info4 integer , dimension(:),allocatable, save :: cap1_info integer , dimension(:),allocatable, save :: cap2_info integer , dimension(:),allocatable, save :: cap3_info integer , dimension(:),allocatable, save :: cap4_info integer , dimension(:),allocatable, save :: cap1_ipvt integer , dimension(:),allocatable, save :: cap2_ipvt integer , dimension(:),allocatable, save :: cap3_ipvt integer , dimension(:),allocatable, save :: cap4_ipvt integer , dimension(:,:),allocatable, save :: ipvt1 integer , dimension(:,:),allocatable, save :: ipvt2 integer , dimension(:,:),allocatable, save :: ipvt3 integer , dimension(:,:),allocatable, save :: ipvt4 double precision, dimension(:,:),allocatable, save :: ainit double precision, dimension(:,:),allocatable, save :: binit ! ! Variables that are only used once ! double precision f1(1:N**2),f2(1:N**2),f3(1:N**2) double precision f4(1:N**2) double precision B1(1:N,1:5) double precision B2(1:N,1:5) double precision C1(1:N,1:3) double precision C2(1:N,1:3) double precision Diag1(1:N) double precision Diag2(1:N) double precision V1(1:2) double precision V2(1:2) double precision Tut1(1:N),Tut2(1:N) double precision konst1(1:N),konst2(1:N) double precision heloes(1:4*N**2) double precision nyshoyre(0:2*N+3,0:2*N+3) double precision sh(0:2*N-1,0:2*N-1) double precision la(0:2*N+3,0:2*N+3) double precision DB1(1:N),DB2(1:N) double precision norm1,norm2,norm3,norm4 ! ! ! Variables to get the capacitance matrix on vector form. ! ! double precision res_Cap1(1:2*N,1:2*N) double precision res_Cap2(1:2*N,1:2*N) double precision res_Cap3(1:2*N,1:2*N) double precision res_Cap4(1:2*N,1:2*N) ! ! A counter is introduced to assure that the necessary ! variables only are allocated once. ! count = 0 count_max = 0 count_rate = 0 call system_clock(count,count_rate,count_max) ocount = count ! ! ! flag = 0 for a first right hand side, if not ! one has a subsequent right hand side. ! if (flag .EQ. 0 .AND. oldn .NE. 0) then ! ! We have been called before, but with a different N, must ! deallocate all arrays before allocating new ones. ! deallocate (Id1) deallocate (Id2) deallocate (Id1trans) deallocate (Id2trans) deallocate (Ja1) deallocate (Ja2) deallocate (Fa1) deallocate (Fa2) deallocate (Sum1) deallocate (Sum2) deallocate (Stor1Sep) deallocate (Stor2Sep) deallocate (Stor3Sep) deallocate (Stor4Sep) deallocate (Cap1) deallocate (Cap2) deallocate (Cap3) deallocate (Cap4) deallocate (De1) deallocate (De2) deallocate (De3) deallocate (De4) deallocate (d) deallocate (r1) deallocate (r2) deallocate (r3) deallocate (s2) deallocate (s3) deallocate (cap1_vekt) deallocate (cap2_vekt) deallocate (cap3_vekt) deallocate (cap4_vekt) deallocate (info1) deallocate (info2) deallocate (info3) deallocate (info4) deallocate (cap1_info) deallocate (cap2_info) deallocate (cap3_info) deallocate (cap4_info) deallocate (cap1_ipvt) deallocate (cap2_ipvt) deallocate (cap3_ipvt) deallocate (cap4_ipvt) deallocate (ipvt1) deallocate (ipvt2) deallocate (ipvt3) deallocate (ipvt4) deallocate (ainit) deallocate (binit) ! end if ! ! dealloacation done ! ! oldn = n ! if (flag .EQ. 0) then ! ! Allocation of the necessary variables, only to be done for ! the first right hand side. ! allocate (Id1(1:N,1:N)) allocate (Id2(1:N,1:N)) allocate (Id1trans(1:N,1:N)) allocate (Id2trans(1:N,1:N)) allocate (Ja1(1:N,1:2)) allocate (Ja2(1:N,1:2)) allocate (Fa1(1:2,1:N)) allocate (Fa2(1:2,1:N)) allocate (Sum1(1:N,1:N)) allocate (Sum2(1:N,1:N)) allocate (Stor1Sep(1:7,1:N,1:N)) allocate (Stor2Sep(1:7,1:N,1:N)) allocate (Stor3Sep(1:7,1:N,1:N)) allocate (Stor4Sep(1:7,1:N,1:N)) allocate (Cap1(1:2*N,1:2*N)) allocate (Cap2(1:2*N,1:2*N)) allocate (Cap3(1:2*N,1:2*N)) allocate (Cap4(1:2*N,1:2*N)) allocate (De1(1:N**2)) allocate (De2(1:N**2)) allocate (De3(1:N**2)) allocate (De4(1:N**2)) allocate (d(0:2*N-1)) allocate (r1(0:2*N-1)) allocate (r2(0:2*N-1)) allocate (r3(0:2*N-1)) allocate (s2(0:2*N-1)) allocate (s3(0:2*N-1)) allocate (cap1_vekt(1:(2*N*(2*N+1)/2))) allocate (cap2_vekt(1:(2*N*(2*N+1)/2))) allocate (cap3_vekt(1:(2*N*(2*N+1)/2))) allocate (cap4_vekt(1:(2*N*(2*N+1)/2))) allocate (info1(1:N)) allocate (info2(1:N)) allocate (info3(1:N)) allocate (info4(1:N)) allocate (cap1_info(1:1)) allocate (cap2_info(1:1)) allocate (cap3_info(1:1)) allocate (cap4_info(1:1)) allocate (cap1_ipvt(1:2*N)) allocate (cap2_ipvt(1:2*N)) allocate (cap3_ipvt(1:2*N)) allocate (cap4_ipvt(1:2*N)) allocate (ipvt1(1:N,1:N)) allocate (ipvt2(1:N,1:N)) allocate (ipvt3(1:N,1:N)) allocate (ipvt4(1:N,1:N)) allocate (ainit(0:2*N+3,0:2*N+3)) allocate (binit(0:2*N+3,0:2*N+3)) ! ! ! The spectral matrices B1,B2,C1,C2 are constructed along with the ! matrices V1 and V2 such that V1V1^T = B1-C1^2 and V2V2^T = B2-C2**2. ! ! odde = 1 call leg_disk(B1,C1,V1,odde,N) odde = 0 call leg_disk(B2,C2,V2,odde,N) ! ! ! The eigenvalue problems Id1^T*C1*Id1=Diag1 and ! Id2^T*C2*Id2=Diag2 are solved and Ja1 = Id1^T*V1 and ! Ja2 = Id2^T*V2 are found. ! ! call leg_eig(C1,V1,Diag1,Id1,Ja1,N) call leg_eig(C2,V2,Diag2,Id2,Ja2,N) ! ! ! ! Compute the first part of P, E_i^T G_i^(-1) E_i , i=1,2 ! ! call leg_pfirst(Ja1,Diag1,alf,bet,Sum1,N,x2) call leg_pfirst(Ja2,Diag2,alf,bet,Sum2,N,x2) ! Id1trans = transpose(Id1) Id2trans = transpose(Id2) ! Fa1 = transpose(Ja1) Fa2 = transpose(Ja2) ! do i = 1,N DB1(i) = Diag1(i)**2 DB2(i) = Diag2(i)**2 end do ! ! ! Factorization of the block diagonal matrix hatZ from (4.16). ! ! call leg_sepfakt(Stor1Sep,B1,C1,DB1,Diag1,alf,bet,info1,ipvt1,N,x1,x2) call leg_sepfakt(Stor2Sep,B2,C2,DB1,Diag1,alf,bet,info2,ipvt2,N,x1,x2) call leg_sepfakt(Stor3Sep,B1,C1,DB2,Diag2,alf,bet,info3,ipvt3,N,x1,x2) call leg_sepfakt(Stor4Sep,B2,C2,DB2,Diag2,alf,bet,info4,ipvt4,N,x1,x2) ! ! Tut1 = 0.5d0*(1.d0+alf*Diag1) ! Tut2 = 0.5d0*(1.d0+alf*Diag2) ! Tut1 = 0.5d0*(x2**2*1.d0+x2*alf*Diag1) Tut2 = 0.5d0*(x2**2*1.d0+x2*alf*Diag2) ! ! ! Construction of the matrices D^(k) from (4.37) ! ! konst1 = 0.5d0*(x1**2*1.d0+bet*DB1+alf*x1*Diag1) konst2 = 0.5d0*(x1**2*1.d0+bet*DB2+alf*x1*Diag2) ! ii = 0 do i = 1,N**2,N ii = ii+1 ! konst1 = 0.5*( 1.d0+bet*DB1(ii)+alf*Diag1(ii) ) ! konst2 = 0.5*( 1.d0+bet*DB2(ii)+alf*Diag2(ii) ) De1(i:i+N-1) = 0.5d0/( konst1(ii)*DB1+DB1(ii)*Tut1+x1*x2*Diag1(ii)*Diag1 ) De2(i:i+N-1) = 0.5d0/( konst1(ii)*DB2+DB1(ii)*Tut2+x1*x2*Diag1(ii)*Diag2 ) De3(i:i+N-1) = 0.5d0/( konst2(ii)*DB1+DB2(ii)*Tut1+x1*x2*Diag2(ii)*Diag1 ) De4(i:i+N-1) = 0.5d0/( konst2(ii)*DB2+DB2(ii)*Tut2+x1*x2*Diag2(ii)*Diag2 ) end do ! ! ! The capacitance matrices are computed. ! ! Cap1 = leg_capmatr(Ja1,Ja1,Diag1,De1,Sum1,alf,bet,N,x1) Cap2 = leg_capmatr(Ja2,Ja1,Diag1,De2,Sum2,alf,bet,N,x1) Cap3 = leg_capmatr(Ja1,Ja2,Diag2,De3,Sum1,alf,bet,N,x1) Cap4 = leg_capmatr(Ja2,Ja2,Diag2,De4,Sum2,alf,bet,N,x1) ! res_Cap1 = Cap1 res_Cap2 = Cap2 res_Cap3 = Cap3 res_Cap4 = Cap4 ! ! ! The capacitance matrices are factorized by a symmetric ! factorization. ! ! call dpotrf('U', 2*N, Cap1, 2*N, info) cap1_info = info if (info.ne.0) then write(6,*)'Error in dpotrf (1),info=',info ! ! ! The capacitance matrix is indefinite and is ! transformed to vektor form in order to use dsptrf. ! ! do i = 1,2*N do ii = i,2*N cap1_vekt(i+(ii-1)*ii/2) = res_Cap1(i,ii) end do end do ! call dsptrf('U',2*N,cap1_vekt,cap1_ipvt,info) if (info.ne.0) write(6,*)'Error in dsptrf, info =',info end if ! call dpotrf('U', 2*N, Cap2, 2*N, info) cap2_info = info if(info.ne.0) then write(6,*)'Error in dpotrf (2),info=',info ! ! ! The capactance matrix is indefinite and is transformed ! to vektor form in order to use dsptrf. ! ! do i = 1,2*N do ii = i,2*N cap2_vekt(i+(ii-1)*ii/2) = res_Cap2(i,ii) end do end do ! call dsptrf('U',2*N,cap2_vekt,cap2_ipvt,info) if (info.ne.0) write(6,*)'Error in dsptrf, info =',info end if ! call dpotrf('U', 2*N, Cap3, 2*N, info) cap3_info = info if(info.ne.0) then write(6,*)'Error in dpotrf (3),info=',info ! ! ! The capacitance matrix is indefinte and is transformed ! to vektor form in order to use dsptrf. ! ! do i = 1,2*N do ii = i,2*N cap3_vekt(i+(ii-1)*ii/2) = res_Cap3(i,ii) end do end do ! call dsptrf('U',2*N,cap3_vekt,cap3_ipvt,info) if (info.ne.0) write(6,*)'Error in dsptrf, info =',info end if ! call dpotrf('U', 2*N, Cap4, 2*N, info) cap4_info = info if(info.ne.0) then write(6,*)'Error in dpotrf (4),info=',info ! ! ! The capacitance matrix is indefinite and is transformed to ! vektor form in order to use dsptrf. ! ! do i = 1,2*N do ii = i,2*N cap4_vekt(i+(ii-1)*ii/2) = res_Cap4(i,ii) end do end do ! call dsptrf('U',2*N,cap4_vekt,cap4_ipvt,info) if (info.ne.0) write(6,*)'Error in dsptrf, info =',info end if ! ! ! Initialization for the Legendre transformation ! ! call legcoef(ainit,binit,xjac,2*N+3) ! ! ! The coefficients needed to go from Legendre basis ! to Shen's special basis are computed. ! ! do i = 0,2*N-1 k=i d(i) = 1.d0/sqrt(2*(2*k+3.d0)**2*(2*k+5.d0)) r1(i) = 2.d0/(2*k+1.d0) r2(i) = -4.d0/(2*k+7.d0) r3(i) = 2*(2*k+3.d0)/((2*k+7.d0)*(2*k+9.d0)) s2(i) = -2*(2*k+5.d0)/(2*k+7.d0) s3(i) = (2*k+3)/(2*k+7.d0) end do ! ! ! End of calculations that only need to be done for ! a first right hand side. ! ! endif ! ! Preprocessing time ! call system_clock(count,count_rate,count_max) tid1 = (count-ocount)/(1.d0 * count_rate) ocount = count ! ! ! Transformation from physical to spectral space. ! ! call phtosp2(hoyre,nyshoyre,binit,2*N+3) call system_clock(count,count_rate,count_max) tid3=(count-ocount)/(1.d0 * count_rate) ocount = count ! ! ! Transformation from Legendre basis to Shen's special basis. ! ! do i=0,2*N-1 ! 112 N^2 flops do ii = 0,2*N-1 hoyre(i,ii) = d(i)*d(ii)*(nyshoyre(i,ii)*r1(i)*r1(ii)+ & nyshoyre(i,ii+2)*r1(i)*r2(ii)+ & nyshoyre(i,ii+4)*r1(i)*r3(ii)+ & nyshoyre(i+2,ii)*r2(i)*r1(ii)+ & nyshoyre(i+2,ii+2)*r2(i)*r2(ii)+ & nyshoyre(i+2,ii+4)*r2(i)*r3(ii)+ & nyshoyre(i+4,ii)*r3(i)*r1(ii)+ & nyshoyre(i+4,ii+2)*r3(i)*r2(ii)+ & nyshoyre(i+4,ii+4)*r3(i)*r3(ii)) end do end do ! ! ! The right-hand side is split into four separate parts. ! ! ! ! ii = 0 do l = 1,N do i = 1,N ii = ii+1 f1(ii) = hoyre(2*(i-1),2*(l-1)) f2(ii) = hoyre(2*i-1,2*(l-1)) f3(ii) = hoyre(2*(i-1),2*l-1) f4(ii) = hoyre(2*i-1,2*l-1) end do end do ! ! 100 format(F20.10) ! ! ! ! ! Norms of f1,...,f4 are computed to test for zero right hand side. ! ! norm1 = 0.d0 norm2 = 0.d0 norm3 = 0.d0 norm4 = 0.d0 ! do i = 1,N**2 norm1 = norm1+f1(i)**2 norm2 = norm2+f2(i)**2 norm3 = norm3+f3(i)**2 norm4 = norm4+f4(i)**2 end do ! norm1 = sqrt(norm1) norm2 = sqrt(norm2) norm3 = sqrt(norm3) norm4 = sqrt(norm4) ! ! ! The four systems M_1 u_1 = f_1 , ..., M_4 u_4 = f_4 from ! (4.22)-(4.25) are solved using leg_hele. ! ! if (norm1.gt.rhs_norm) then call leg_hele(Stor1Sep,Id1,Id1trans,Id1,Id1trans, & Cap1,cap1_vekt,Ja1,Fa1,f1,info1,ipvt1,cap1_info, & cap1_ipvt,N) else f1 = 0.d0 end if ! ! if (norm2.gt.rhs_norm) then call leg_hele(Stor2Sep,Id1,Id1trans,Id2,Id2trans, & Cap2,cap2_vekt,Ja1,Fa1,f2,info2,ipvt2,cap2_info, & cap2_ipvt,N) ! else ! f2 = 0.d0 ! end if ! ! ! if (norm3.gt.rhs_norm) then call leg_hele(Stor3Sep,Id2,Id2trans,Id1,Id1trans, & Cap3,cap3_vekt,Ja2,Fa2,f3,info3,ipvt3,cap3_info, & cap3_ipvt,N) ! else ! f3 = 0.d0 ! end if ! ! if (norm4.gt.rhs_norm) then call leg_hele(Stor4Sep,Id2,Id2trans,Id2,Id2trans, & Cap4,cap4_vekt,Ja2,Fa2,f4,info4,ipvt4,cap4_info, & cap4_ipvt,N) ! else ! f4 = 0.d0 ! end if ! ! ! The four solutions are merged together. ! ! ii = 0 do i = 1,N do l = 1,N ii = ii+1 heloes(4*N*(i-1)+2*l-1) = f1(ii) heloes(4*N*(i-1)+2*l) = f2(ii) heloes(4*N*i-2*N+2*l-1) = f3(ii) heloes(4*N*i-2*N+2*l) = f4(ii) end do end do ! ! ii = 0 do i = 1,4*N**2,2*N sh(:,ii) = heloes(i:i+2*N-1) ii = ii+1 end do ! ! ! The solution is transformed from Shen's basis back to ! Legendre basis. ! ! call basistrans(la,sh,d,s2,s3,N) ! coeff = la ! print '(F4.2)',4.77 ! print '(F8.3)',la ! ! Time inside Legendre solver ! call system_clock(count,count_rate,count_max) tid2 = (count-ocount)/(1.d0 * count_rate) ocount = count ! ! ! Transformation of the solution from spectral space to physical space. ! ! call sptoph2(hoyre,la,ainit,2*N+3) call system_clock(count,count_rate,count_max) tid3 = tid3 + (count-ocount)/(1.d0 * count_rate) ! tid4 = tid1+tid2+tid3 ! end subroutine leg