! ! CODE to model the contribution of pericytes and macrophages in ! the initiation of angiogenesis. Ref: Levine, Sleeman & Nilsen-Hamilton ! Quadratic elements in space and backward Euler in time ! ! UNKNOWNS (in order) ! u1 (nuk=1) - chemotactic agent ! v1 (nuk=2) - angiogenic factor ! c1 (nuk=3) - protease ! f1 (nuk=4) - fibronectin ! a1 (nuk=5) - angiostatic agent ! pi1 (nuk=6) - protease inhibitor ! sig1 (nuk=7) - pericyte cell density ! rmac1 (nuk=8) - macrophage cell density ! eta1 (nuk=9) - endothelial cell density ! ! *********************************************************************** ! ! implicit double precision (a-h,o-z) ! Grid parameters integer,parameter:: nx=81 ! number of nodes in x direction (not including midpoints) ! ! parameters for 1-d calculation integer,parameter:: ldn1d=161 ! for dimensioning (number of nodes) integer,parameter:: lde1d=80 ! for dimensioning (number of elements) integer,parameter:: lda1d=3 ! for dimensioning (leading dimension of matrices) integer,parameter:: nlocn1d=3 ! number of local nodes in quadratic element integer,parameter:: nq1d=3 ! 3 point Gauss quadrature rule used integer,parameter:: neqn1d=9 ! number of equations integer,parameter:: maxiter_cap=30 ! maximum number of iterations in capillary double precision, parameter:: tol_cap=0.5d-8 ! tolerance for convergence in capillary ! ! arrays for 1-d calculation - geometry double precision xc1d(1:ldn1d) ! x-coordinate of node double precision qx1d(1:lde1d,1:nq1d) ! quadrature points per element double precision qw1d(1:lde1d,1:nq1d) ! quadrature weight per element double precision area1d(1:lde1d) ! area of each element integer indx1d(1:ldn1d) ! unknown number integer node1d(1:lde1d,1:nlocn1d) ! local node information ! ! arrays for 1-d calculation - solution double precision amass1d(1:lda1d,1:ldn1d) ! mass matrix double precision alap1d(1:lda1d,1:ldn1d) ! mass + dt*stiffness matrix double precision u1old(1:ldn1d) ! in capillary at previous timestep double precision u1(1:ldn1d) ! in capillary double precision v1old(1:ldn1d) ! angiogenic factor in capillary at previous timestep double precision v1(1:ldn1d) ! angiogenic factor in capillary double precision c1old(1:ldn1d) ! proteolytic enzyme in capillary at previous timestep double precision c1(1:ldn1d) ! proteolytic enzyme in capillary double precision f1old(1:ldn1d) ! fibronectin in capillary at previous timestep double precision f1(1:ldn1d) ! fibronectin in capillary double precision a1old(1:ldn1d) ! in capillary at previous timestep double precision a1(1:ldn1d) ! in capillary double precision pi1old(1:ldn1d) ! in capillary at previous timestep double precision pi1(1:ldn1d) ! in capillary double precision sig1old(1:ldn1d) ! in capillary at previous timestep double precision sig1(1:ldn1d) ! in capillary double precision rmac1old(1:ldn1d) ! in capillary at previous timestep double precision rmac1(1:ldn1d) ! in capillary double precision eta1old(1:ldn1d) ! endothelial cell density in capillary at previous timestep double precision eta1(1:ldn1d) ! endothelial cell density in capillary ! ! parameters for time stepping integer,parameter:: init=0 ! =0 set initial conditions at t=0; =1 read from file integer,parameter:: maxit_time=10000 ! maximum number of time steps double precision,parameter:: time_final=1d0 ! final time double precision,parameter:: dtinit=0.0005d0 ! initial timestep double precision,parameter:: dt_reducefactor=0.5d0 ! reduction factor for decreasing dt double precision,parameter:: dt_incfactor=1.0d0 ! factor for decreasing dt ! ! parameter for outputting data integer,parameter:: nplot=5 ! write out solution every nplot !timesteps ! include 'perimac.data' ! ! open files for writing ! unit 5 - writes off miscellaneous data open(unit=5,file='output.m') open(unit=2,file='plotdata.m') ncall=0 ! *********************************************************************** ! GEOMETRY ROUTINES FOR CAPILLARY & ECM ! *********************************************************************** call geom1d(xc1d,qx1d,qw1d,area1d,node1d,indx1d, & nel1d,nunk1d, np1d,nq1d,ib1d,nlocn1d, & nx,lde1d ) ! ! ! *********************************************************************** ! ASSEMBLE AND FACTOR MASS AND STIFFNESS COEFFICIENT MATRICES ! *********************************************************************** ! call assmass1d(amass1d,xc1d,qx1d,qw1d,node1d, & indx1d,nlocn1d,nunk1d,nq1d,ib1d,nel1d,lde1d,lda1d ) ! dt=dtinit ! ! ! *********************************************************************** ! SET INITIAL CONDITIONS or READ STARTING VALUES FROM FILE ! *********************************************************************** ! ! set initial guess for solution of equations (init=0) ! or read starting values from file (init=1) if(init == 0) then time=0.0d0 call set_init1d(u1old,v1old,c1old,f1old,a1old,pi1old, & sig1old,rmac1old,eta1old, & xc1d,indx1d,nunk1d) ! ! set initial guess for unknowns u1(1:nunk1d)=u1old(1:nunk1d) v1(1:nunk1d)=v1old(1:nunk1d) c1(1:nunk1d)=c1old(1:nunk1d) f1(1:nunk1d)=f1old(1:nunk1d) a1(1:nunk1d)=a1old(1:nunk1d) pi1(1:nunk1d)=pi1old(1:nunk1d) sig1(1:nunk1d)=sig1old(1:nunk1d) rmac1(1:nunk1d)=rmac1old(1:nunk1d) eta1(1:nunk1d)=eta1old(1:nunk1d) end if ! ! *********************************************************************** ! ITERATION LOOP AT EACH TIMESTEP ! *********************************************************************** do niter=1,maxit_time ! 101 continue time=time+dt write(*,1021) time write(5,1021) time ! ! ! *********************************************************************** ! step 2. solve equations in capillary ! *********************************************************************** ! iconv_cap=0 ! convergence flag: <=0 converges, =1 failed ! =-1 converged in < 3 iterations call solvecap(iconv_cap,time,dt,alap1d,amass1d, & u1old,u1,v1old,v1,c1old,c1,f1old,f1,a1old,a1, & pi1old,pi1,sig1old,sig1,rmac1old,rmac1,eta1old,eta1, & xc1d,qx1d,qw1d,node1d,indx1d,nlocn1d,nunk1d,nq1d,ib1d,& nel1d,tol_cap,maxiter_cap,neqn1d,lde1d,lda1d ) ! ! if iteration failed to converge in capillary, then reduce timestep ! and begin again if( iconv_cap > 0 ) then go to 401 else go to 501 end if ! ! *********************************************************************** ! iteration failed to converge so reduce timestep ! *********************************************************************** ! 401 continue time=time-dt dt=dt_reducefactor * dt write(*,1031) dt u1(1:nunk1d)=u1old(1:nunk1d) v1(1:nunk1d)=v1old(1:nunk1d) c1(1:nunk1d)=c1old(1:nunk1d) f1(1:nunk1d)=f1old(1:nunk1d) a1(1:nunk1d)=a1old(1:nunk1d) pi1(1:nunk1d)=pi1old(1:nunk1d) sig1(1:nunk1d)=sig1old(1:nunk1d) rmac1(1:nunk1d)=rmac1old(1:nunk1d) eta1(1:nunk1d)=eta1old(1:nunk1d) go to 101 ! ! *********************************************************************** ! Get ready for next time step ! ! check to see if timestep can be increased ! check to see if reached final time ! *********************************************************************** ! 501 continue u1old(1:nunk1d)=u1(1:nunk1d) v1old(1:nunk1d)=v1(1:nunk1d) c1old(1:nunk1d)=c1(1:nunk1d) f1old(1:nunk1d)=f1(1:nunk1d) a1old(1:nunk1d)=a1(1:nunk1d) pi1old(1:nunk1d)=pi1(1:nunk1d) sig1old(1:nunk1d)=sig1(1:nunk1d) rmac1old(1:nunk1d)=rmac1(1:nunk1d) eta1old(1:nunk1d)=eta1(1:nunk1d) ! if (iconv_cap < 0 ) then dt=dt*dt_incfactor end if ! if (time >= time_final) then call output(time, ncall, & u1, v1, c1,f1, a1,pi1,sig1,rmac1,eta1,xc1d,indx1d, nx,np1d, neqn1d ) stop end if ! ! *********************************************************************** ! write off data ! *********************************************************************** if ( mod(niter,nplot) == 0 ) call output( & time, ncall, u1, v1, c1,f1, a1,pi1,sig1,rmac1,eta1, & xc1d,indx1d, nx,np1d, neqn1d ) ! end do ! *********************************************************************** stop ! 1021 format(/' time = ',f12.5) 1031 format(/,' reducing the timestep to ', e12.5) end ! ! ! *********************************************************************** ! subroutine asslap1d(nuk,dt,a,xc,qx,qw,node,indx,nlocn, & nunk,nq,ib,nel,lde,lda ) ! ! *********************************************************************** ! Sets up the 1-D coefficient matrix for u'v'+dt*uv ! ! INPUT ! xc, qx, qw,area - coordinate and quadrature info from geometry ! indx,node,nlocn,nunk,nq,ib,nel - info from geometry ! ! OUTPUT ! a - the factored coefficient matrix ! ! ! *********************************************************************** ! implicit double precision (a-h,o-z) dimension a(lda,1), xc(1), qx(lde,1), qw(lde,1), & indx(*),node(lde,1) ! ! include 'perimac.data' ! ! zero matrix a(1:lda,1:nunk)=0.0d0 ! ! set diffusivity constant if (nuk == 7 ) dif = cap_dif_sig if (nuk == 8 ) dif = cap_dif_mac if (nuk == 9 ) dif = cap_dif_eta ! ! assemble matrix and right hand side do it=1,nel ! do iq=1,nq x=qx(it,iq) ar=qw(it,iq) ! do in=1,nlocn ip=node(it,in) i=indx(ip ) ! if(i > 0) then call qbf1d(x,it,in,bb,bx,xc,node,lde) do inn=1,nlocn ipp=node(it,inn) j=indx(ipp) ! ! matrix is symmetric so just form lower part if(j > 0) then if(j <= i ) then call qbf1d(x,it,inn,bbb,bbx,xc,node,lde) ij=i-j+1 aij= bb*bbb+dt*dif*bx*bbx a(ij,j)=a(ij,j)+aij*ar end if end if ! end do ! end if ! end do ! end do ! end do ! ! factor matrix and store info=0 call DPBTRF( 'L',nunk,ib,a,lda,info) if(info /= 0 ) then write(*,1001) info write(5,1001) info stop end if ! return ! 1001 format(' info not equal to zero in factoring laplacian ',i10) end ! ! ! *********************************************************************** ! *********************************************************************** subroutine assmass1d(a,xc,qx,qw,node,indx, & nlocn,nunk,nq,ib,nel,lde,lda ) ! *********************************************************************** ! Sets up the 1-D mass matrix and factors ! ! INPUT ! xc,qx,qw,area - coordinate and quadrature info from geometry ! indx,node,nlocn,nunk,nq,ib,nel - info from geometry ! ! OUTPUT ! a - the factored mass matrix ! ! *********************************************************************** ! implicit double precision (a-h,o-z) dimension a(lda,1), xc(1), qx(lde,1),qw(lde,1), & indx(1),node(lde,1) ! ! zero matrix a(1:lda,1:nunk)=0.0d0 ! ! assemble mass matrix element by element ! do it=1,nel ! do iq=1,nq x=qx(it,iq) ar=qw(it,iq) ! do in=1,nlocn ip=node(it,in) i=indx(ip ) if(i > 0) then call qbf1d(x,it,in,bb,bx,xc,node,lde) ! do inn=1,nlocn ipp=node(it,inn) j=indx(ipp) ! ! matrix is symmetric so just form lower part if(j <= i) then call qbf1d(x,it,inn,bbb,bbx,xc,node,lde) aij=bb*bbb ij=i-j+1 a(ij,j)=a(ij,j)+aij*ar end if ! end do ! end if ! end do end do end do ! ! ! factor mass matrix and store info=0 CALL DPBTRF( 'L',nunk,ib,a,lda,info) if(info.ne.0) then write(*,1010) info write(5,1010) info stop endif ! return ! 1010 format(' matrix factor failed in assmass and info=',i10) end ! ! *********************************************************************** subroutine geom1d(xc,qx,qw,area, node,indx, & nel,nunk, np,nq,ib, nlocn, nx, lde ) ! *********************************************************************** ! Sets up geometric information for a rectangular domain ! ! INPUT: ! nlocn - number of local nodes for quadratic element (3) ! nq - number of quadrature points for chosen rule ! xl,xr - left and right endpoints of domain ! lde - leading dimensions of element arrays ! ! OUTPUT: ! xc - x coordinates of nodes ! qx - x coordinates of quadrature points ! qw - quadrature weights ! area - area of each element ! node - array which associates to each local node of element the ! corresponding global node ! indx - array associating unknown number to global node ! nel - number of elements ! nunk - number of unknowns ! np - number of nodes ! ib - half bandwidth of coefficient matrix ! ! *********************************************************************** ! implicit double precision (a-h,o-z) dimension xc(1), area(1),qx(lde,1),qw(lde,1), & node(lde,1),indx(*),xq(5),wq(5) ! ! mapping functions for variable grid grdx(x)=( (xr-xl)/2.0d0 ) / sqrt((xr-xl)/2.0d0)*sqrt(x) ! include 'perimac.data' ! ! set up grid information and unknown counter array nxm=nx-1 hx=(xr-xl)/dfloat(nxm) nx2=nx+nxm iu=0 ip=0 it1=0 xx=xl-hx/2.d0 ! do ic=1,nx2 xx=xx+hx/2.d0 ix=mod(ic,2) ip=ip+1 ! if(ic.le.nx) then xc(ip)=grdx(xx) else if (ic > nx.and.ic < nx2) then xc(ip)=xr-grdx(xr-xx) else if(ic == nx2) then xc(ip)=xr endif ! ! Don't set up node value if at midpoint or last node if(ix /= 1 .or. ic == nx2 ) go to 120 it1=it1+1 node(it1,1)=ip node(it1,2)=ip+1 node(it1,3)=ip+ 2 120 continue ! ! Neumann boundary conditions iu=iu+1 indx(ip)=iu 125 continue end do ! nunk=iu np=ip nel=it1 write(*,1001) nunk,np,nel write(5,1001) nunk,np,nel ! ! reset midpoint node to midpoint do it=1,nel ip1=node(it,1) ip2=node(it,2) ip3=node(it,3) x1=xc(ip1) x3=xc(ip3) xc(ip2)=(x1+x3)/2.0d0 end do ! write(5,1005) (i,xc(i),i=1,np) write(5,1006) (i,indx(i) ,i=1,np) write(5,1007) (i,node(i,1),node(i,2),node(i,3),i=1,nel) ! ! set up quadrature information needed for assembly call gauss3pt(xq,wq) ! do it=1,nel ip1=node(it,1) ip3=node(it,3) x1=xc(ip1) x3=xc(ip3) area(it)=abs(x3-x1) x1=xc( node(it,1) ) ! qx(it,1:nq)=x1+area(it)/2.d0*(xq(1:nq)+1.d0) qw(it,1:nq)=0.5d0*area(it)*wq(1:nq) ! end do ! write(5,1011) (i,qx(i,1),qx(i,2),qx(i,3), i=1,nel ) write(5,1012) (i,qw(i,1),qw(i,2),qw(i,3), i=1,nel ) write(5,1013) (i,area(i),i=1,nel) ! ! ! determine bandwidth of coefficient matrix for nu unknowns ib=0 do it=1,nel ! do in=1,nlocn ip=node(it,in) i=indx(ip) if(i > 0) then ! do inn=1,nlocn ipp=node(it,inn) j=indx(ipp) ij=j-i if(ij.gt.ib) ib=ij end do end if ! end do ! end do ! ibt=ib+ib+1 write(*,1009) ibt write(5,1009) ibt ! return ! 1001 format(/,' Geometry for Capillary',/ & ,' number of unknowns=',i5,/' number of nodes=',i5,/,& ' number of elements=',i5) 1005 format(' xc', i10, f12.6) 1006 format(' indx',6i10 ) 1007 format(' node array',4i10) 1009 format(/' total bandwidth for 1D problem',i5) 1011 format(' quad pt',i10,3f12.5) 1012 format(' quad weight', i10,3f12.5) 1013 format(' area of element', i10,f12.5) end ! ! *********************************************************************** ! subroutine resid1d(residn,time,dt, & u1old,u1,v1old,v1,c1old,c1,f1old,f1,a1old,a1,pi1old,pi1, & sig1old,sig1,rmac1old,rmac1,eta1old,eta1, & xc,qx,qw,node,indx, & nlocn,nunk,nq,nel,neqn1d,lde ) ! ! *********************************************************************** ! ! Sets up the residual vector in capillary ! ! *********************************************************************** ! implicit double precision (a-h,o-z) ! dimension u1(*),u1old(*),v1(*),v1old(*),c1(*),c1old(*),f1(*),f1old(*), & a1(*),a1old(*),pi1(*),pi1old(*),sig1(*),sig1old(*), & rmac1(*),rmac1old(*),eta1(*),eta1old(*), & xc(1),qx(lde,1),qw(lde,1),node(lde,1),indx( *), & re(2000) ! include 'perimac.data' ! ! zero residual array ndimre=neqn1d*nunk re(1:ndimre)=0.0d0 ! ! assemble residual ! do it=1,nel ! do iq=1,nq ar=qw(it,iq) x=qx(it,iq) ! ! calculate value of solution from previous time step at quad point call evalpt1d(x,uhold,uholdx,u1old,xc,indx,node,it,nlocn,lde ) call evalpt1d(x,vhold,vholdx,v1old,xc,indx,node,it,nlocn,lde ) call evalpt1d(x,chold,choldx,c1old,xc,indx,node,it,nlocn,lde ) call evalpt1d(x,fhold,fholdx,f1old,xc,indx,node,it,nlocn,lde ) call evalpt1d(x,ahold,aholdx,a1old,xc,indx,node,it,nlocn,lde ) call evalpt1d(x,pihold,piholdx,pi1old,xc,indx,node,it,nlocn,lde ) call evalpt1d(x,sighold,sigholdx,sig1old,xc,indx,node,it,nlocn,lde ) call evalpt1d(x,rmachold,rmacholdx,rmac1old,xc,indx,node,it,nlocn,lde ) call evalpt1d(x,etahold,etaholdx,eta1old,xc,indx,node,it,nlocn,lde ) ! ! evaluate current solution at quad point call evalpt1d(x,uh,uhx,u1,xc,indx,node,it,nlocn,lde ) call evalpt1d(x,vh,vhx,v1,xc,indx,node,it,nlocn,lde ) call evalpt1d(x,ch,chx,c1,xc,indx,node,it,nlocn,lde ) call evalpt1d(x,fh,fhx,f1,xc,indx,node,it,nlocn,lde ) call evalpt1d(x,ah,ahx,a1,xc,indx,node,it,nlocn,lde ) call evalpt1d(x,pih,pihx,pi1,xc,indx,node,it,nlocn,lde ) call evalpt1d(x,sigh,sighx,sig1,xc,indx,node,it,nlocn,lde ) call evalpt1d(x,rmach,rmachx,rmac1,xc,indx,node,it,nlocn,lde ) call evalpt1d(x,etah,etahx,eta1,xc,indx,node,it,nlocn,lde ) cah=ch / (1.0d0 + cap_nue * pih + cap_nu4 * fh) ! ! assemble residual do in=1,nlocn ip=node(it,in) i=indx(ip) ! if (i > 0) then call qbf1d(x,it,in,bb,bx,xc,node,lde) ! do nuk=1,neqn1d ! if (nuk == 1) then term_time=bb*(uh-uhold)/dt term_de=bb*cap_lambda2*uh*rmach/(1.0d0+cap_nu2*uh) term_rhs=bb*rhs1d_source(nuk,time,x,xc,indx,node,it,nlocn,lde) else if (nuk == 2) then term_time=bb*(vh-vhold)/dt term_de=-bb*cap_lambda2*uh*rmach/(1.0d0+cap_nu2*uh) & + bb*cap_lambda1*vh*etah/(1.0d0+cap_nu1*vh) term_rhs=bb*rhs1d_source(nuk,time,x,xc,indx,node,it,nlocn,lde) else if (nuk == 3) then term_time=bb*(ch-chold)/dt term_de=bb*( -cap_lambda1*vh*etah/(1.0d0+cap_nu1*vh) + cap_mu*ch ) term_rhs=bb*rhs1d_source(nuk,time,x,xc,indx,node,it,nlocn,lde) else if (nuk == 4) then term_time=bb*(fh-fhold)/dt term_de=bb*( -cap_beta*(cap_fzero-fh)*fh*etah+cap_lambda4*cah*fh/ & (1.0d0+cap_nu4*fh) ) term_rhs=bb*rhs1d_source(nuk,time,x,xc,indx,node,it,nlocn,lde) else if (nuk == 5) then term_time=bb*(ah-ahold)/dt term_de=cap_lambda3*ah*etah/(1.0d0+cap_nu3*ah) *bb term_rhs=bb*rhs1d_source(nuk,time,x,xc,indx,node,it,nlocn,lde) else if (nuk == 6) then term_time=bb*(pih-pihold)/dt term_de=-bb*(cap_lambda3*ah*etah/(1.0d0+cap_nu3*ah) ) term_rhs=bb*rhs1d_source(nuk,time,x,xc,indx,node,it,nlocn,lde) else if (nuk == 7) then term_time=bb*(sigh-sighold)/dt call trans_prob_sig(fh,fhx, tauxovertau) term_de= cap_dif_sig*( sighx -sigh*tauxovertau )*bx term_rhs=bb*rhs1d_source(nuk,time,x,xc,indx,node,it,nlocn,lde) else if (nuk == 8) then term_time=bb*(rmach-rmachold)/dt call trans_prob_mac(uh,uhx,tauxovertau) term_de= cap_dif_mac*( rmachx -rmach*tauxovertau )*bx term_rhs=bb*rhs1d_source(nuk,time,x,xc,indx,node,it,nlocn,lde) else if (nuk == 9) then term_time=bb*(etah-etahold)/dt call trans_prob_eta(fh,fhx,ch,chx,pih,pihx,tauxovertau) term_de= cap_dif_eta*( etahx -etah*tauxovertau )*bx term_rhs=bb*rhs1d_source(nuk,time,x,xc,indx,node,it,nlocn,lde) end if ! k=i+(nuk-1)*nunk re(k)=re(k)+(term_time+term_de-term_rhs)*ar end do ! end if ! end do ! end do ! end do ! ! ! calculate norm of residual residn=sum( ( re(1:ndimre) )**2) residn=sqrt(residn)/dfloat(ndimre) write(*,1001) residn ! return ! 1001 format(' norm of residual ',e12.5) end ! ! ! *********************************************************************** subroutine solvecap(iconv,time,dt,alap,am, & u1old,u1,v1old,v1,c1old,c1,f1old,f1,a1old,a1, & pi1old,pi1,sig1old,sig1,rmac1old,rmac1,eta1old,eta1, & xc,qx,qw,node,indx, & nlocn,nunk,nq,ib,nel,tol,maxiter,neqn1d,lde,lda) ! *********************************************************************** ! ! solves the problem in the capillary given the values in ECM ! ! INPUT ! a,am factored coefficient matrices ! f - work vector ! vold,...eold - values of solution at previous timestep ! v1...e1 - initial guesses at new timestep(values at old timestep) ! xc,area,node,indx,nlocn,nunk,nq,ib,nel - info from geometry ! maxiter- maximum number of iterations ! toln - tolerance for convergence ! time,dt - timestepping info ! lde,lda,ldn - dimensioning info ! ! OUTPUT ! v1...e1 - converged value of solution at this timestep ! ! *********************************************************************** ! implicit double precision (a-h,o-z) ! dimension alap(lda,1),am(lda,1), & u1old(*),u1(*),v1old(*),v1(*),c1old(*),c1(*),f1old(*),f1(*), & a1old(*),a1(*),pi1old(*),pi1(*),sig1old(*),sig1(*), & rmac1old(*),rmac1(*),eta1old(*),eta1(*), & xc(1),qx(lde,1),qw(lde,1),node(lde,1),indx( *), & f(200) ! ! ! *********************************************************************** ! ITERATION LOOP FOR CONVERGENCE IN CAPILLARY ! *********************************************************************** ! residn_old=1.0d5 do niter=1,maxiter ! ! *********************************************************************** ! solve for u - chemotactic agent equation - nuk= 1 ! *********************************************************************** ! nuk=1 call solsys1d(nuk,time,dt,am,f,u1old,u1,v1,c1,f1,a1,pi1, & sig1,rmac1,eta1,xc,qx,qw, & node,indx,nlocn,nunk,nq,ib,nel,lde,lda) ! u1(1:nunk)=f(1:nunk) ! write(*,*) 'nuk=', nuk do i=1,nunk ! write(*,*) i,u1(i) end do ! ! *********************************************************************** ! solve for angiogenic factor equation - nuk= 2 ! *********************************************************************** ! nuk=2 call solsys1d(nuk,time,dt,am,f,v1old,u1,v1,c1,f1,a1,pi1, & sig1,rmac1,eta1,xc,qx,qw, & node,indx,nlocn,nunk,nq,ib,nel,lde,lda) ! v1(1:nunk)=f(1:nunk) ! write(*,*) 'nuk=', nuk do i=1,nunk ! write(*,*) i,v1(i) end do ! ! *********************************************************************** ! solve for protolytic enzyme equation - nuk = 3 ! *********************************************************************** ! nuk=3 call solsys1d(nuk,time,dt,am,f,c1old,u1,v1,c1,f1,a1,pi1, & sig1,rmac1,eta1,xc,qx,qw, & node,indx,nlocn,nunk,nq,ib,nel,lde,lda) ! c1(1:nunk)=f(1:nunk) ! write(*,*) 'nuk=', nuk do i=1,nunk ! write(*,*) i,c1(i) end do ! ! ! *********************************************************************** ! solve for fibronectin equation - nuk = 4 ! *********************************************************************** ! nuk=4 call solsys1d(nuk,time,dt,am,f,f1old,u1,v1,c1,f1,a1,pi1, & sig1,rmac1,eta1,xc,qx,qw, & node,indx,nlocn,nunk,nq,ib,nel,lde,lda) ! f1(1:nunk)=f(1:nunk) ! write(*,*) 'nuk=', nuk do i=1,nunk ! write(*,*) i,f1(i) end do ! ! ! *********************************************************************** ! solve for angiostatic agent equation - nuk = 5 ! *********************************************************************** ! nuk=5 call solsys1d(nuk,time,dt,am,f,a1old,u1,v1,c1,f1,a1,pi1, & sig1,rmac1,eta1,xc,qx,qw, & node,indx,nlocn,nunk,nq,ib,nel,lde,lda) ! a1(1:nunk)=f(1:nunk) ! write(*,*) 'nuk=', nuk do i=1,nunk ! write(*,*) i,a1(i) end do ! ! ! *********************************************************************** ! solve for protease inhibitor equation - nuk = 6 ! *********************************************************************** ! nuk=6 call solsys1d(nuk,time,dt,am,f,pi1old,u1,v1,c1,f1,a1,pi1, & sig1,rmac1,eta1,xc,qx,qw, & node,indx,nlocn,nunk,nq,ib,nel,lde,lda) ! pi1(1:nunk)=f(1:nunk) ! write(*,*) 'nuk=', nuk do i=1,nunk ! write(*,*) i,pi1(i) end do ! ! ! *********************************************************************** ! solve for PC density equation - nuk = 7 ! *********************************************************************** ! nuk=7 call asslap1d(nuk,dt,alap,xc,qx,qw,node, & indx,nlocn,nunk,nq,ib,nel,lde,lda) ! call solsys1d(nuk,time,dt,alap,f,sig1old,u1,v1,c1,f1,a1,pi1, & sig1,rmac1,eta1,xc,qx,qw, & node,indx,nlocn,nunk,nq,ib,nel,lde,lda) ! sig1(1:nunk)=f(1:nunk) ! write(*,*) 'nuk=', nuk do i=1,nunk ! write(*,*) i,sig1(i) end do ! ! ! *********************************************************************** ! solve for macrophages equation - nuk = 8 ! *********************************************************************** ! nuk=8 call asslap1d(nuk,dt,alap,xc,qx,qw,node, & indx,nlocn,nunk,nq,ib,nel,lde,lda) ! call solsys1d(nuk,time,dt,alap,f,rmac1old,u1,v1,c1,f1,a1,pi1, & sig1,rmac1,eta1,xc,qx,qw, & node,indx,nlocn,nunk,nq,ib,nel,lde,lda) ! rmac1(1:nunk)=f(1:nunk) ! write(*,*) 'nuk=', nuk do i=1,nunk ! write(*,*) i,rmac1(i) end do ! ! ! *********************************************************************** ! solve for endothelial cell density equation - nuk = 9 ! *********************************************************************** ! nuk=9 call asslap1d(nuk,dt,alap,xc,qx,qw,node, & indx,nlocn,nunk,nq,ib,nel,lde,lda) ! call solsys1d(nuk,time,dt,alap,f,eta1old,u1,v1,c1,f1,a1,pi1, & sig1,rmac1,eta1,xc,qx,qw, & node,indx,nlocn,nunk,nq,ib,nel,lde,lda) ! eta1(1:nunk)=f(1:nunk) ! write(*,*) 'nuk=', nuk do i=1,nunk ! write(*,*) i,eta1(i) end do ! ! ! *********************************************************************** ! calculate residual for system and check for convergence ! *********************************************************************** ! call resid1d(residn,time,dt, & u1old,u1, v1old,v1,c1old,c1,f1old,f1, & a1old,a1,pi1old,pi1,sig1old,sig1,rmac1old,rmac1,eta1old,eta1, & xc,qx,qw,node,indx,nlocn,nunk,nq,nel,neqn1d,lde ) ! if( residn < tol ) then write(*,1001) niter if ( niter < 3) iconv=-1 return end if ! ! if(niter > 1 ) then ! if (residn >= residn_old) then ! if residual not decreasing, stop ! iconv=1 ! return ! end if residn_old=residn end if ! end do ! ! iteration failed to converge write(*,1002) iconv=1 ! return ! 1001 format(' solution converged in capillary in ',i5,' iterations') 1002 format(' solution failed to converge in capillary ') end ! ! ! *********************************************************************** ! subroutine solsys1d(nuk,time,dt,a,f,fold,u1,v1,c1,f1,a1,pi1, & sig1,rmac1,eta1, & xc,qx,qw,node,indx,nlocn,nunk,nq,ib,nel,lde,lda ) ! ! *********************************************************************** ! ! Sets up the right-hand side vector using ! piecewise quadratics and solves system given L L^t of given matrix ! ! INPUT ! nuk - which unknown we are solving for ! time,dt - timestepping info ! amat - correct coefficient matrix (am or alap) ! fold - value of current unknown at previous timestep (for rhs) ! v1...e1 - current values of variables ! xc,yc,node,indx,nlocn,nunk,nq,ib,nel - geometry variables ! lde,lda,ldn - dimensioning variables ! xc2d,yc2d,indx2d,node2d,lde2d - info to use value of V in ECM as ! source term in equation ! ! OUTPUT ! f - solution of unknown # nuk ! ! *********************************************************************** ! implicit double precision (a-h,o-z) dimension a(lda,1),fold(1),f(1), & u1(*),v1(1),c1(1),f1(1),a1(*),pi1(*),sig1(*),rmac1(*),eta1(1),& xc(1),qx(lde,1),qw(lde,1),indx(1),node(lde,1) ! ! ! zero rhs f(1:nunk)=0.0d0 ! ! assemble right hand side do it=1,nel ! do iq=1,nq ar=qw(it,iq) x=qx(it,iq) ! ! evaluate source term on right hand side of equation rhs_source= rhs1d_source(nuk,time,x,xc,indx,node,it,nlocn,lde ) ! ! evaluate solution at previous time step for time derivative call evalpt1d(x,soln,solnx,fold,xc,indx,node,it,nlocn,lde ) ! do in=1,nlocn ip=node(it,in) i=indx(ip) ! if(i > 0) then call qbf1d(x,it,in,bb,bx,xc,node,lde) ! ! evaluate right hand side of equation consisting of lagged terms rhs_lag= rhs1d_lag(nuk,time,x,bb,bx, & u1,v1,c1,f1,a1,pi1,sig1,rmac1,eta1, & xc,indx,node,it,nlocn,lde ) f(i)=f(i)+ar*( bb*(soln +dt*rhs_source) + dt*rhs_lag ) end if ! end do ! end do ! end do ! ! solve system nrhs=1 info=0 call dpbtrs('L',nunk,ib,nrhs,a,lda,f,nunk,info) if(info /= 0) then write(*,1001) info stop end if return 1001 format(' info in banded Cholesky solver is',i10) end ! ! *********************************************************************** ! *********************************************************************** ! ! USER SUPPLIED ROUTINES ! ! rhs1d_source - gives source term on right hand side of equations ! rhs1d_lag - gives lagged terms on right hand side of equations ! solninit1d - initial conditions for 1-D problem ! transprob1d - transition probability function for eta equation ! ! *********************************************************************** ! ! *********************************************************************** function rhs1d_source(nuk,time,x,xc,indx,node,it,nlocn,lde ) ! *********************************************************************** ! ! Gives source term on right-hand side of DE ! ! *********************************************************************** ! implicit double precision (a-h,o-z) ! dimension xc(*),indx(*),node(lde,*) ! include 'perimac.data' ! pi=4.d0*atan(1.0d0) ! ! ! u1 equation if (nuk == 1) then term=cap_kappa*( 1.d0-cos( 2.d0*pi*x/(xr-xl) ) ) rhs1d_source= cap_uzero* (term )**cap_mexp * exp(-cap_theta*time) return end if ! v1 equation if (nuk == 2) then term = cap_kappa*(1.d0-cos(2.d0*pi*x/(xr-xl) ) ) rhs1d_source= cap_vzero* (term )**cap_mexp * exp(-cap_theta*time) return end if ! ! c1 equation if (nuk == 3) then rhs1d_source= 0.0d0 return end if ! ! f1 equation if (nuk == 4) then rhs1d_source= 0.0d0 return end if ! ! ! a1 equation if (nuk == 5) then rhs1d_source = cap_ar return end if ! ! pi1 equation if (nuk == 6) then rhs1d_source= 0.0d0 return end if ! ! sig1 equation if (nuk == 7) then rhs1d_source= 0.0d0 return end if ! ! mac1 equation if (nuk == 8) then rhs1d_source= 0.0d0 return end if ! eta1 equation if (nuk == 9) then rhs1d_source=0.0d0 return end if return end ! ! ! *********************************************************************** function rhs1d_lag(nuk,time,x,bb,bx, & u1,v1,c1,f1,a1,pi1,sig1,rmac1,eta1,xc,indx,node,it,nlocn,lde ) ! *********************************************************************** ! ! Gives lagged terms on the right-hand side of DE ! terms are multiplied by appropriate basis function ! ! *********************************************************************** ! implicit double precision (a-h,o-z) ! dimension u1(*),v1(*),c1(*),f1(*),a1(*),sig1(*),rmac1(*),eta1(*), & xc(*),indx(*),node(lde,*) ! include 'perimac.data' ! ! u1 equation if (nuk == 1) then call evalpt1d(x,uh,uhx,u1,xc,indx,node,it,nlocn,lde ) call evalpt1d(x,rmach,rmachx,rmac1,xc,indx,node,it,nlocn,lde ) rhs1d_lag=( -cap_lambda2*uh*rmach/(1.0d0+cap_nu2*uh) ) * bb return end if ! ! v1 equation if (nuk == 2) then call evalpt1d(x,vh,vhx,v1,xc,indx,node,it,nlocn,lde ) call evalpt1d(x,etah,etahx,eta1,xc,indx,node,it,nlocn,lde ) call evalpt1d(x,uh,uhx,u1,xc,indx,node,it,nlocn,lde ) call evalpt1d(x,rmach,rmachx,rmac1,xc,indx,node,it,nlocn,lde ) rhs1d_lag=(cap_lambda2*uh*rmach/(1.0d0+cap_nu2*uh) & -cap_lambda1*vh*etah/(1.0d0+cap_nu1*vh) ) * bb return end if ! ! c1 equation if (nuk == 3) then call evalpt1d(x,vh,vhx,v1,xc,indx,node,it,nlocn,lde ) call evalpt1d(x,etah,etahx,eta1,xc,indx,node,it,nlocn,lde ) call evalpt1d(x,ch,chx,c1,xc,indx,node,it,nlocn,lde ) rhs1d_lag=( cap_lambda1*vh*etah/(1.0d0+cap_nu1*vh)-cap_mu*ch ) *bb return end if ! ! f1 equation if (nuk == 4) then call evalpt1d(x,fh,fhx,f1,xc,indx,node,it,nlocn,lde ) call evalpt1d(x,etah,etahx,eta1,xc,indx,node,it,nlocn,lde ) call evalpt1d(x,ch,chx,c1,xc,indx,node,it,nlocn,lde ) call evalpt1d(x,pih,pihx,pi1,xc,indx,node,it,nlocn,lde ) cah=ch/(1.0d0 + cap_nue * pih + cap_nu4 * fh) rhs1d_lag=( cap_beta*(cap_fzero-fh)*fh*etah & -cap_lambda4*cah*fh/(1.0d0+cap_nu4*fh) ) *bb return end if ! ! a1 equation if (nuk == 5) then call evalpt1d(x,ah,ahx,a1,xc,indx,node,it,nlocn,lde ) call evalpt1d(x,etah,etahx,eta1,xc,indx,node,it,nlocn,lde ) rhs1d_lag=( -cap_lambda3*ah*etah/(1.0d0+cap_nu3*ah) ) * bb return end if ! pi1 equation if (nuk == 6) then call evalpt1d(x,ah,ahx,a1,xc,indx,node,it,nlocn,lde ) call evalpt1d(x,etah,etahx,eta1,xc,indx,node,it,nlocn,lde ) rhs1d_lag=( cap_lambda3*ah*etah/(1.0d0+cap_nu3*ah) ) * bb return end if ! sig1 equation if (nuk == 7) then call evalpt1d(x,fh,fhx,f1,xc,indx,node,it,nlocn,lde ) call evalpt1d(x,sigh,sighx,sig1,xc,indx,node,it,nlocn,lde ) call trans_prob_sig(fh,fhx, tauxovertau) rhs1d_lag=cap_dif_sig*( sigh*tauxovertau ) *bx return end if ! rmac1 equation if (nuk == 8) then call evalpt1d(x,uh,uhx,u1,xc,indx,node,it,nlocn,lde ) call evalpt1d(x,rmach,rmachx,rmac1,xc,indx,node,it,nlocn,lde ) call trans_prob_mac(uh,uhx,tauxovertau) rhs1d_lag=cap_dif_mac*(rmach*tauxovertau ) *bx return end if ! eta1 equation if (nuk == 9) then call evalpt1d(x,fh,fhx,f1,xc,indx,node,it,nlocn,lde ) call evalpt1d(x,etah,etahx,eta1,xc,indx,node,it,nlocn,lde ) call evalpt1d(x,ch,chx,c1,xc,indx,node,it,nlocn,lde ) call evalpt1d(x,pih,pihx,pi1,xc,indx,node,it,nlocn,lde ) call trans_prob_eta(fh,fhx,ch,chx,pih,pihx,tauxovertau) rhs1d_lag=cap_dif_eta*( etah*tauxovertau ) *bx return end if ! return end ! ! ! ! *********************************************************************** subroutine trans_prob_eta(fh,fhx,ch,chx,pih,pihx,tauxovertau) ! *********************************************************************** ! ! Given f,pin,c and a the constants alpha1,alpha2, beta1, beta2 & gamma1, ! gamma2 and nu (rnut) this computes ! tau and its derivatives ! ! *********************************************************************** implicit double precision (a-h,o-z) ! include 'perimac.data' ! term= 1.0d0 + cap_nue * pih + cap_nu4 * fh cah=ch/term cahx=chx/term-ch*(cap_nue*pihx + cap_nu4*fhx) / (term**2) delta_alpha=cap_alpha2 - cap_alpha1 delta_beta=cap_beta2 - cap_beta1 termc=cap_gamma1 *( cahx*(delta_alpha) ) / ((cah+cap_alpha1 )*(cah+cap_alpha2 ) ) termf=cap_gamma2 *( fhx*(delta_beta) )/( (fh+cap_beta1 )*(fh+cap_beta2 ) ) tauxovertau=termc+termf return end ! ! ! *********************************************************************** subroutine trans_prob_sig(fh,fhx, tauxovertau) ! *********************************************************************** ! ! Given f,pin,c and a the constants alpha1,alpha2, beta1, beta2 & gamma1, ! gamma2 and nu (rnut) this computes ! tau and its derivatives ! ! *********************************************************************** implicit double precision (a-h,o-z) ! include 'perimac.data' ! delta_alpha=cap_alpha4 - cap_alpha3 prod= (fh+cap_alpha3) * (fh+cap_alpha4) tauxovertau= cap_gamma3 * fhx * delta_alpha / prod return end ! ! ! *********************************************************************** subroutine trans_prob_mac(uh,uhx,tauxovertau) ! *********************************************************************** ! ! Given f,pin,c and a the constants alpha1,alpha2, beta1, beta2 & gamma1, ! gamma2 and nu (rnut) this computes ! tau and its derivatives ! ! *********************************************************************** implicit double precision (a-h,o-z) ! include 'perimac.data' ! delta_beta=cap_beta4 - cap_beta3 prod= (uh+cap_beta3) * (uh+cap_beta4) tauxovertau= cap_gamma4 * uhx * delta_beta / prod return end ! ! *********************************************************************** subroutine set_init1d(u1,v1,c1,f1,a1,pi1,sig1,rmac1,eta1,xc,indx,np) ! *********************************************************************** ! ! Gives the initial condition for time dependent problem ! Must be supplied by user ! ! ! *********************************************************************** implicit double precision (a-h,o-z) dimension u1(*),v1(*),c1(*),f1(*),a1(*),pi1(*), & sig1(*),rmac1(*),eta1(*),xc(*),indx(*) ! ! include 'perimac.data' ! pi=4.*atan(1.) ! do ip=1,np i=indx(ip) if (i > 0 ) then u1(i)=0.0d0 v1(i)=0.0d0 c1(i)=0.0d0 f1(i)=1.0d0 a1(i)=0.0d0 pi1(i)=0.0d0 sig1(i)=1.0d0 rmac1(i)=1.0d0 eta1(i)=1.0d0 end if end do ! return end ! *********************************************************************** ! *********************************************************************** ! ! UTILITY ROUTINES ! ! gauss3pt - 3 point Gauss rule in 1-D ! qpt1d - quadratic basis functions in 1-D ! ! *********************************************************************** ! ! *********************************************************************** ! subroutine evalpt1d(x,soln,solnx,f,xc,indx,node,it,nlocn,lde ) ! ! *********************************************************************** ! evaluates the 1-D solution in capillary & its derivatives at a point x ! given its nodal values in f ! ! INPUT: ! x - x point where solution is to be evaluated ! f - vector of nodal values (f) ! xc - x-coordinates of nodes ! indx & node arrays from geometry ! nlocn - the number of nodes in each element ! OUTPUT: soln,solnx - the value & derivative at point x ! ! *********************************************************************** ! implicit double precision (a-h,o-z) ! dimension f(*),xc(*),indx(*),node(lde,*) ! soln=0.d0 solnx=0.d0 ! do in=1,nlocn ip=node(it,in) call qbf1d(x,it,in,bb,bx,xc,node,lde) i=indx(ip) ! if( i > 0 ) then soln =soln +bb*f(i) solnx=solnx+bx*f(i) end if ! end do ! return end ! ! *********************************************************************** ! subroutine gauss3pt(xq,wq) ! ! *********************************************************************** implicit double precision(a-h,o-z) ! dimension xq(1), wq(1) ! xq(1)=-.7745966692d0 xq(2)=0.0d0 xq(3)=-xq(1) wq(1)=0.5555555556d0 wq(2)=0.8888888889d0 wq(3)=wq(1) ! return end ! ! ! *********************************************************************** ! subroutine gauss5pt(xq,wq) ! ! *********************************************************************** implicit double precision(a-h,o-z) ! dimension xq(1), wq(1) ! xq(1)=.9061798459d0 xq(2)=.5384693101d0 xq(3)=0.0d0 xq(4)=-xq(2) xq(5)=-xq(1) wq(1)=0.2369268850d0 wq(2)=0.4786286705d0 wq(3)=0.5688888889d0 wq(4)=wq(2) wq(5)=wq(1) ! return end ! ! ! *********************************************************************** ! *********************************************************************** subroutine output ( time, ncall, u1, v1, c1,f1, a1,pi1,sig1,rmac1,eta1, & xc, indx, nx, np, neqn1d ) ! ! *********************************************************************** ! ! Outputs solution to file for matlab plotting ! ! *********************************************************************** ! implicit double precision (a-h,o-z) ! dimension u1(*), v1(*),c1(*),f1(*),& a1(*),pi1(*),sig1(*),rmac1(*),eta1(*), & xc(*), indx(*) ! include 'perimac.data' ! ncall=ncall+1 npx=nx+nx-1 ! ! write grid information if first call to output ! if(ncall.eq.1) then write(2,1010) write(2,1020) ( xc(i),i=1,np-1 ) ! x-coordinates of vertices write(2,1011) xc(np) end if ! ! write out time ! if ( (mod(ncall, 2) == 0).OR.(ncall < 2) ) then !jmj write(2,1040) ncall, time ! ! ******************************************************************* ! write out solution in capillary at nodes of elements ! ******************************************************************* ! do nwrite=1,neqn1d ! do ip=1,np i=indx(ip) ! write u1 if( nwrite == 1 ) then if(ip == 1) write(2,2010) ncall if(ip /= np) write(2,1060) u1(i) if(ip == np) write(2,1061) u1(i) ! write v1 else if( nwrite == 2 ) then if(ip == 1) write(2,1070) ncall if(ip /= np) write(2,1060) v1(i) if(ip == np) write(2,1061) v1(i) ! write c1 else if(nwrite == 3) then if(ip == 1) write(2,1080) ncall if(ip /= np) write(2,1060) c1(i) if(ip == np) write(2,1061) c1(i) ! write f1 else if(nwrite == 4) then if(ip == 1) write(2,1090) ncall if(ip /= np) write(2,1060) f1(i) if(ip == np) write(2,1061) f1(i) ! write a1 else if( nwrite == 5 ) then if(ip == 1) write(2,2020) ncall if(ip /= np) write(2,1060) a1(i) if(ip == np) write(2,1061) a1(i) ! write pi1 else if( nwrite == 6 ) then if(ip == 1) write(2,2030) ncall if(ip /= np) write(2,1060) pi1(i) if(ip == np) write(2,1061) pi1(i) ! write sig1 else if( nwrite == 7 ) then if(ip == 1) write(2,2040) ncall if(ip /= np) write(2,1060) sig1(i) if(ip == np) write(2,1061) sig1(i) ! write rmac1 else if( nwrite == 8 ) then if(ip == 1) write(2,2050) ncall if(ip /= np) write(2,1060) rmac1(i) if(ip == np) write(2,1061) rmac1(i) ! write e1 else if(nwrite == 9) then if(ip == 1) write(2,1050) ncall if(ip /= np) write(2,1060) eta1(i) if(ip == np) write(2,1061) eta1(i) endif ! end do ! end loop over nodes ! end do ! end loop over unknowns ! ! end if !jmj ! 1010 format("x =[") 1011 format(e12.5,"];" ) 1020 format(e12.5,";",e12.5,";",e12.5,";",e12.5,";", & e12.5,";",e12.5,";" ) 1040 format("t(",i5,")=",f12.5,";") 1050 format("eta1(:,",i5,")=[") 1060 format(e12.5,";") 1061 format(e12.5,"];") 1070 format("v1(:,",i5,")=[") 1080 format("c1(:,",i5,")=[") 1090 format("f1(:,",i5,")=[") 2010 format("u1(:,",i5,")=[") 2020 format("a1(:,",i5,")=[") 2030 format("pi1(:,",i5,")=[") 2040 format("sig1(:,",i5,")=[") 2050 format("rmac1(:,",i5,")=[") end ! ! *********************************************************************** subroutine qbf1d(x,it,in,bb,bx,xc,node,lde) ! *********************************************************************** ! Quadratic basis functions in 1-D ! INPUT: ! x - point, it - triangle, in - node basis function centered at ! xc - x coordinates of grid node - triangle info lde - dimensioning ! OUTPUT: ! bb, bx - basis function and derivative at point ! *********************************************************************** ! implicit double precision (a-h,o-z) dimension node(lde,1),xc(1) ! in1=in in2=mod(in,3)+1 in3=mod(in+1,3)+1 i1=node(it,in1) i2=node(it,in2) i3=node(it,in3) den=(xc(i3)-xc(i1))*(xc(i2)-xc(i1)) ! if(den.eq.0) then write(*,1001) i1,i2,i3 stop end if ! bb=(xc(i3)-x)*(xc(i2)-x)/den bx=-(xc(i2)-x)/den-(xc(i3)-x)/den ! return 1001 format(' den in qbf is zero ',3i10 ) end ! ! ! *********************************************************************** SUBROUTINE DPBSV( UPLO, N, KD, NRHS, AB, LDAB, B, LDB, INFO ) ! ! -- LAPACK driver routine (version 1.0) -- ! Univ. of Tennessee, Univ. of California Berkeley, NAG Ltd., ! Courant Institute, Argonne National Lab, and Rice University ! February 29, 1992 ! ! .. Scalar Arguments .. CHARACTER UPLO INTEGER INFO, KD, LDAB, LDB, N, NRHS ! .. ! .. Array Arguments .. DOUBLE PRECISION AB( LDAB, * ), B( LDB, * ) ! .. ! ! Purpose ! ======= ! ! DPBSV computes the solution to a real system of linear equations ! A * X = B, ! where A is an N by N symmetric positive definite band matrix and X ! and B are N by NRHS matrices. ! ! The Cholesky decomposition is used to factor A as ! A = U'* U , if UPLO = 'U', or ! A = L * L', if UPLO = 'L', ! where U is an upper triangular matrix, L is a lower triangular ! matrix, and ' indicates transpose. The factored form of A ! is then used to solve the system of equations A * X = B. ! ! Arguments ! ========= ! ! UPLO (input) CHARACTER*1 ! Specifies whether the upper or lower triangular part of the ! symmetric matrix A is stored: ! = 'U': Upper triangular ! = 'L': Lower triangular ! ! N (input) INTEGER ! The number of linear equations, i.e., the order of the ! matrix A. N >= 0. ! ! KD (input) INTEGER ! The number of super-diagonals of the matrix A if UPLO = 'U', ! or the number of sub-diagonals if UPLO = 'L'. KD >= 0. ! ! NRHS (input) INTEGER ! The number of right hand sides, i.e., the number of columns ! of the matrix B. NRHS >= 0. ! ! AB (input/output) DOUBLE PRECISION array, dimension (LDAB,N) ! On entry, the upper or lower triangle of the symmetric band ! matrix A, stored in the first KD+1 rows of the array. The ! j-th column of A is stored in the j-th column of the array AB ! as follows: ! if UPLO = 'U', AB(KD+1+i-j,j) = A(i,j) for max(1,j-KD)<=i<=j; ! if UPLO = 'L', AB(1+i-j,j) = A(i,j) for j<=i<=min(N,j+KD). ! See below for further details. ! ! On exit, if INFO = 0, the triangular factor U or L from the ! Cholesky factorization A = U'*U or A = L*L' of the band ! matrix A, in the same storage format as A. ! ! LDAB (input) INTEGER ! The leading dimension of the array AB. LDAB >= KD+1. ! ! B (input/output) DOUBLE PRECISION array, dimension (LDB,NRHS) ! On entry, the N by NRHS matrix of right hand side vectors B ! for the system of equations A*X = B. ! On exit, if INFO = 0, the N by NRHS matrix of solution ! vectors X. ! ! LDB (input) INTEGER ! The leading dimension of the array B. LDB >= max(1,N). ! ! INFO (output) INTEGER ! = 0: successful exit ! < 0: if INFO = -k, the k-th argument had an illegal value ! > 0: if INFO = k, the leading minor of order k of A is not ! positive definite, so the factorization could not be ! completed, and the solution has not been computed. ! ! Further Details ! =============== ! ! The band storage scheme is illustrated by the following example, when ! N = 6, KD = 2, and UPLO = 'U': ! ! On entry: On exit: ! ! * * a13 a24 a35 a46 * * u13 u24 u35 u46 ! * a12 a23 a34 a45 a56 * u12 u23 u34 u45 u56 ! a11 a22 a33 a44 a55 a66 u11 u22 u33 u44 u55 u66 ! ! Similarly, if UPLO = 'L' the format of A is as follows: ! ! On entry: On exit: ! ! a11 a22 a33 a44 a55 a66 l11 l22 l33 l44 l55 l66 ! a21 a32 a43 a54 a65 * l21 l32 l43 l54 l65 * ! a31 a42 a53 a64 * * l31 l42 l53 l64 * * ! ! Array elements marked * are not used by the routine. ! ! ===================================================================== ! ! .. External Functions .. LOGICAL LSAME EXTERNAL LSAME ! .. ! .. External Subroutines .. EXTERNAL DPBTRF, dpbtrs, XERBLA ! .. ! .. Intrinsic Functions .. INTRINSIC MAX ! .. ! .. Executable Statements .. ! ! Test the input parameters. ! INFO = 0 IF( .NOT.LSAME( UPLO, 'U' ) .AND. .NOT.LSAME( UPLO, 'L' ) ) THEN INFO = -1 ELSE IF( N.LT.0 ) THEN INFO = -2 ELSE IF( KD.LT.0 ) THEN INFO = -3 ELSE IF( NRHS.LT.0 ) THEN INFO = -4 ELSE IF( LDAB.LT.KD+1 ) THEN INFO = -6 ELSE IF( LDB.LT.MAX( 1, N ) ) THEN INFO = -8 END IF IF( INFO.NE.0 ) THEN CALL XERBLA( 'DPBSV ', -INFO ) RETURN END IF ! ! Compute the Cholesky factorization A = U'*U or A = L*L'. ! CALL DPBTRF( UPLO, N, KD, AB, LDAB, INFO ) IF( INFO.EQ.0 ) THEN ! ! Solve the system A*X = B, overwriting B with X. ! CALL dpbtrs( UPLO, N, KD, NRHS, AB, LDAB, B, LDB, INFO ) ! END IF RETURN ! ! End of DPBSV ! END SUBROUTINE DPBTRF( UPLO, N, KD, AB, LDAB, INFO ) ! ! -- LAPACK routine (version 1.0) -- ! Univ. of Tennessee, Univ. of California Berkeley, NAG Ltd., ! Courant Institute, Argonne National Lab, and Rice University ! February 29, 1992 ! ! .. Scalar Arguments .. CHARACTER UPLO INTEGER INFO, KD, LDAB, N ! .. ! .. Array Arguments .. DOUBLE PRECISION AB( LDAB, * ) ! .. ! ! Purpose ! ======= ! ! DPBTRF computes the Cholesky factorization of a real symmetric ! positive definite band matrix A. ! ! The factorization has the form ! A = U' * U , if UPLO = 'U', or ! A = L * L', if UPLO = 'L', ! where U is an upper triangular matrix and L is lower triangular. ! ! Arguments ! ========= ! ! UPLO (input) CHARACTER*1 ! Specifies whether the upper or lower triangular part of the ! symmetric matrix A is stored: ! = 'U': Upper triangular ! = 'L': Lower triangular ! ! N (input) INTEGER ! The order of the matrix A. N >= 0. ! ! KD (input) INTEGER ! The number of super-diagonals of the matrix A if UPLO = 'U', ! or the number of sub-diagonals if UPLO = 'L'. KD >= 0. ! ! AB (input/output) DOUBLE PRECISION array, dimension (LDAB,N) ! On entry, the upper or lower triangle of the symmetric band ! matrix A, stored in the first KD+1 rows of the array. The ! j-th column of A is stored in the j-th column of the array AB ! as follows: ! if UPLO = 'U', AB(kd+1+i-j,j) = A(i,j) for max(1,j-kd)<=i<=j; ! if UPLO = 'L', AB(1+i-j,j) = A(i,j) for j<=i<=min(n,j+kd). ! ! On exit, if INFO = 0, the triangular factor U or L from the ! Cholesky factorization A = U'*U or A = L*L' of the band ! matrix A, in the same storage format as A. ! ! LDAB (input) INTEGER ! The leading dimension of the array AB. LDAB >= KD+1. ! ! INFO (output) INTEGER ! = 0: successful exit ! < 0: if INFO = -k, the k-th argument had an illegal value ! > 0: if INFO = k, the leading minor of order k is not ! positive definite, and the factorization could not be ! completed. ! ! Further Details ! =============== ! ! The band storage scheme is illustrated by the following example, when ! N = 6, KD = 2, and UPLO = 'U': ! ! On entry: On exit: ! ! * * a13 a24 a35 a46 * * u13 u24 u35 u46 ! * a12 a23 a34 a45 a56 * u12 u23 u34 u45 u56 ! a11 a22 a33 a44 a55 a66 u11 u22 u33 u44 u55 u66 ! ! Similarly, if UPLO = 'L' the format of A is as follows: ! ! On entry: On exit: ! ! a11 a22 a33 a44 a55 a66 l11 l22 l33 l44 l55 l66 ! a21 a32 a43 a54 a65 * l21 l32 l43 l54 l65 * ! a31 a42 a53 a64 * * l31 l42 l53 l64 * * ! ! Array elements marked * are not used by the routine. ! ! Contributed by ! Peter Mayes and Giuseppe Radicati, IBM ECSEC, Rome, March 23, 1989 ! ! ===================================================================== ! ! .. Parameters .. DOUBLE PRECISION ONE, ZERO PARAMETER ( ONE = 1.0D+0, ZERO = 0.0d0) INTEGER NBMAX, LDWORK PARAMETER ( NBMAX = 32, LDWORK = NBMAX+1 ) ! .. ! .. Local Scalars .. INTEGER I, I2, I3, IB, II, J, JJ, NB ! .. ! .. Local Arrays .. DOUBLE PRECISION WORK( LDWORK, NBMAX ) ! .. ! .. External Functions .. LOGICAL LSAME INTEGER ILAENV EXTERNAL LSAME, ILAENV ! .. ! .. External Subroutines .. EXTERNAL DGEMM, DPBTF2, DPOTF2, DSYRK, DTRSM, XERBLA ! .. ! .. Intrinsic Functions .. INTRINSIC MIN ! .. ! .. Executable Statements .. ! ! Test the input parameters. ! INFO = 0 IF( ( .NOT.LSAME( UPLO, 'U' ) ) .AND. & ( .NOT.LSAME( UPLO, 'L' ) ) ) THEN INFO = -1 ELSE IF( N.LT.0 ) THEN INFO = -2 ELSE IF( KD.LT.0 ) THEN INFO = -3 ELSE IF( LDAB.LT.KD+1 ) THEN INFO = -5 END IF IF( INFO.NE.0 ) THEN CALL XERBLA( 'DPBTRF', -INFO ) RETURN END IF ! ! Quick return if possible ! IF( N.EQ.0 ) & RETURN ! ! Determine the block size for this environment ! NB = ILAENV( 1, 'DPBTRF', UPLO, N, KD, -1, -1 ) ! ! The block size must not exceed the semi-bandwidth KD, and must not ! exceed the limit set by the size of the local array WORK. ! NB = MIN( NB, NBMAX ) ! IF( NB.LE.1 .OR. NB.GT.KD ) THEN ! ! Use unblocked code ! CALL DPBTF2( UPLO, N, KD, AB, LDAB, INFO ) ELSE ! ! Use blocked code ! IF( LSAME( UPLO, 'U' ) ) THEN ! ! Compute the Cholesky factorization of a symmetric band ! matrix, given the upper triangle of the matrix in band ! storage. ! ! Zero the upper triangle of the work array. ! DO 20 J = 1, NB DO 10 I = 1, J - 1 WORK( I, J ) = ZERO 10 CONTINUE 20 CONTINUE ! ! Process the band matrix one diagonal block at a time. ! DO 70 I = 1, N, NB IB = MIN( NB, N-I+1 ) ! ! Factorize the diagonal block ! CALL DPOTF2( UPLO, IB, AB( KD+1, I ), LDAB-1, II ) IF( II.NE.0 ) THEN INFO = I + II - 1 GO TO 150 END IF IF( I+IB.LE.N ) THEN ! ! Update the relevant part of the trailing submatrix. ! If A11 denotes the diagonal block which has just been ! factorized, then we need to update the remaining ! blocks in the diagram: ! ! A11 A12 A13 ! A22 A23 ! A33 ! ! The numbers of rows and columns in the partitioning ! are IB, I2, I3 respectively. The blocks A12, A22 and ! A23 are empty if IB = KD. The upper triangle of A13 ! lies outside the band. ! I2 = MIN( KD-IB, N-I-IB+1 ) I3 = MIN( IB, N-I-KD+1 ) ! IF( I2.GT.0 ) THEN ! ! Update A12 ! CALL DTRSM( 'Left', 'Upper', 'Transpose', & 'Non-unit', IB, I2, ONE, AB( KD+1, I ),& LDAB-1, AB( KD+1-IB, I+IB ), LDAB-1 ) ! ! Update A22 ! CALL DSYRK( 'Upper', 'Transpose', I2, IB, -ONE,& AB( KD+1-IB, I+IB ), LDAB-1, ONE,& AB( KD+1, I+IB ), LDAB-1 ) END IF ! IF( I3.GT.0 ) THEN ! ! Copy the lower triangle of A13 into the work array. ! DO 40 JJ = 1, I3 DO 30 II = JJ, IB WORK( II, JJ ) = AB( II-JJ+1, JJ+I+KD-1 ) 30 CONTINUE 40 CONTINUE ! ! Update A13 (in the work array). ! CALL DTRSM( 'Left', 'Upper', 'Transpose',& 'Non-unit', IB, I3, ONE, AB( KD+1, I ),& LDAB-1, WORK, LDWORK ) ! ! Update A23 ! IF( I2.GT.0 ) & CALL DGEMM( 'Transpose', 'No Transpose', I2, I3,& IB, -ONE, AB( KD+1-IB, I+IB ),& LDAB-1, WORK, LDWORK, ONE,& AB( 1+IB, I+KD ), LDAB-1 ) ! ! Update A33 ! CALL DSYRK( 'Upper', 'Transpose', I3, IB, -ONE,& WORK, LDWORK, ONE, AB( KD+1, I+KD ),& LDAB-1 ) ! ! Copy the lower triangle of A13 back into place. ! DO 60 JJ = 1, I3 DO 50 II = JJ, IB AB( II-JJ+1, JJ+I+KD-1 ) = WORK( II, JJ ) 50 CONTINUE 60 CONTINUE END IF END IF 70 CONTINUE ELSE ! ! Compute the Cholesky factorization of a symmetric band ! matrix, given the lower triangle of the matrix in band ! storage. ! ! Zero the lower triangle of the work array. ! DO 90 J = 1, NB DO 80 I = J + 1, NB WORK( I, J ) = ZERO 80 CONTINUE 90 CONTINUE ! ! Process the band matrix one diagonal block at a time. ! DO 140 I = 1, N, NB IB = MIN( NB, N-I+1 ) ! ! Factorize the diagonal block ! CALL DPOTF2( UPLO, IB, AB( 1, I ), LDAB-1, II ) IF( II.NE.0 ) THEN INFO = I + II - 1 GO TO 150 END IF IF( I+IB.LE.N ) THEN ! ! Update the relevant part of the trailing submatrix. ! If A11 denotes the diagonal block which has just been ! factorized, then we need to update the remaining ! blocks in the diagram: ! ! A11 ! A21 A22 ! A31 A32 A33 ! ! The numbers of rows and columns in the partitioning ! are IB, I2, I3 respectively. The blocks A21, A22 and ! A32 are empty if IB = KD. The lower triangle of A31 ! lies outside the band. ! I2 = MIN( KD-IB, N-I-IB+1 ) I3 = MIN( IB, N-I-KD+1 ) ! IF( I2.GT.0 ) THEN ! ! Update A21 ! CALL DTRSM( 'Right', 'Lower', 'Transpose',& 'Non-unit', I2, IB, ONE, AB( 1, I ),& LDAB-1, AB( 1+IB, I ), LDAB-1 ) ! ! Update A22 ! CALL DSYRK( 'Lower', 'No Transpose', I2, IB, -ONE,& AB( 1+IB, I ), LDAB-1, ONE,& AB( 1, I+IB ), LDAB-1 ) END IF ! IF( I3.GT.0 ) THEN ! ! Copy the upper triangle of A31 into the work array. ! DO 110 JJ = 1, IB DO 100 II = 1, MIN( JJ, I3 ) WORK( II, JJ ) = AB( KD+1-JJ+II, JJ+I-1 ) 100 CONTINUE 110 CONTINUE ! ! Update A31 (in the work array). ! CALL DTRSM( 'Right', 'Lower', 'Transpose',& 'Non-unit', I3, IB, ONE, AB( 1, I ),& LDAB-1, WORK, LDWORK ) ! ! Update A32 ! IF( I2.GT.0 ) & CALL DGEMM( 'No transpose', 'Transpose', I3, I2,& IB, -ONE, WORK, LDWORK,& AB( 1+IB, I ), LDAB-1, ONE,& AB( 1+KD-IB, I+IB ), LDAB-1 ) ! ! Update A33 ! CALL DSYRK( 'Lower', 'No Transpose', I3, IB, -ONE,& WORK, LDWORK, ONE, AB( 1, I+KD ),& LDAB-1 ) ! ! Copy the upper triangle of A31 back into place. ! DO 130 JJ = 1, IB DO 120 II = 1, MIN( JJ, I3 ) AB( KD+1-JJ+II, JJ+I-1 ) = WORK( II, JJ ) 120 CONTINUE 130 CONTINUE END IF END IF 140 CONTINUE END IF END IF RETURN ! 150 CONTINUE RETURN ! ! End of DPBTRF ! END SUBROUTINE dpbtrs( UPLO, N, KD, NRHS, AB, LDAB, B, LDB, INFO ) ! ! -- LAPACK routine (version 1.0) -- ! Univ. of Tennessee, Univ. of California Berkeley, NAG Ltd., ! Courant Institute, Argonne National Lab, and Rice University ! February 29, 1992 ! ! .. Scalar Arguments .. CHARACTER UPLO INTEGER INFO, KD, LDAB, LDB, N, NRHS ! .. ! .. Array Arguments .. DOUBLE PRECISION AB( LDAB, * ), B( LDB, * ) ! .. ! ! Purpose ! ======= ! ! DPBTRS solves a system of linear equations A*X = B with a symmetric ! positive definite band matrix A using the Cholesky factorization ! A = U'*U or A = L*L' computed by DPBTRF. ! ! Arguments ! ========= ! ! UPLO (input) CHARACTER*1 ! Specifies whether the factor stored in AB is upper or lower ! triangular. ! = 'U': Upper triangular ! = 'L': Lower triangular ! ! N (input) INTEGER ! The order of the matrix A. N >= 0. ! ! KD (input) INTEGER ! The number of super-diagonals of the matrix A if UPLO = 'U', ! or the number of sub-diagonals if UPLO = 'L'. KD >= 0. ! ! NRHS (input) INTEGER ! The number of right hand sides, i.e., the number of columns ! of the matrix B. NRHS >= 0. ! ! AB (input) DOUBLE PRECISION array, dimension (LDAB,N) ! The triangular factor U or L from the Cholesky factorization ! A = U'*U or A = L*L' of the band matrix A, stored in the ! first KD+1 rows of the array. The j-th column of U or L is ! stored in the array AB as follows: ! if UPLO ='U', AB(kd+1+i-j,j) = U(i,j) for max(1,j-kd)<=i<=j; ! if UPLO ='L', AB(1+i-j,j) = L(i,j) for j<=i<=min(n,j+kd). ! ! LDAB (input) INTEGER ! The leading dimension of the array AB. LDAB >= KD+1. ! ! B (input/output) DOUBLE PRECISION array, dimension (LDB,NRHS) ! On entry, the right hand side vectors B for the system of ! linear equations. ! On exit, the solution vectors, X. ! ! LDB (input) INTEGER ! The leading dimension of the array B. LDB >= max(1,N). ! ! INFO (output) INTEGER ! = 0: successful exit ! < 0: if INFO = -k, the k-th argument had an illegal value ! ! ===================================================================== ! ! .. Local Scalars .. LOGICAL UPPER INTEGER J ! .. ! .. External Functions .. LOGICAL LSAME EXTERNAL LSAME ! .. ! .. External Subroutines .. EXTERNAL DTBSV, XERBLA ! .. ! .. Intrinsic Functions .. INTRINSIC MAX ! .. ! .. Executable Statements .. ! ! Test the input parameters. ! INFO = 0 UPPER = LSAME( UPLO, 'U' ) IF( .NOT.UPPER .AND. .NOT.LSAME( UPLO, 'L' ) ) THEN INFO = -1 ELSE IF( N.LT.0 ) THEN INFO = -2 ELSE IF( KD.LT.0 ) THEN INFO = -3 ELSE IF( NRHS.LT.0 ) THEN INFO = -4 ELSE IF( LDAB.LT.KD+1 ) THEN INFO = -6 ELSE IF( LDB.LT.MAX( 1, N ) ) THEN INFO = -8 END IF IF( INFO.NE.0 ) THEN CALL XERBLA( 'dpbtrs', -INFO ) RETURN END IF ! ! Quick return if possible ! IF( N.EQ.0 .OR. NRHS.EQ.0 )& RETURN ! IF( UPPER ) THEN ! ! Solve A*X = B where A = U'*U. ! DO 10 J = 1, NRHS ! ! Solve U'*X = B, overwriting B with X. ! CALL DTBSV( 'Upper', 'Transpose', 'Non-unit', N, KD, AB,& LDAB, B( 1, J ), 1 ) ! ! Solve U*X = B, overwriting B with X. ! CALL DTBSV( 'Upper', 'No transpose', 'Non-unit', N, KD, AB,& LDAB, B( 1, J ), 1 ) 10 CONTINUE ELSE ! ! Solve A*X = B where A = L*L'. ! DO 20 J = 1, NRHS ! ! Solve L*X = B, overwriting B with X. ! CALL DTBSV( 'Lower', 'No transpose', 'Non-unit', N, KD, AB,& LDAB, B( 1, J ), 1 ) ! ! Solve L'*X = B, overwriting B with X. ! CALL DTBSV( 'Lower', 'Transpose', 'Non-unit', N, KD, AB,& LDAB, B( 1, J ), 1 ) 20 CONTINUE END IF ! RETURN ! ! End of DPBTRS ! END SUBROUTINE DPBTF2( UPLO, N, KD, AB, LDAB, INFO ) ! ! -- LAPACK routine (version 1.0) -- ! Univ. of Tennessee, Univ. of California Berkeley, NAG Ltd., ! Courant Institute, Argonne National Lab, and Rice University ! February 29, 1992 ! ! .. Scalar Arguments .. CHARACTER UPLO INTEGER INFO, KD, LDAB, N ! .. ! .. Array Arguments .. DOUBLE PRECISION AB( LDAB, * ) ! .. ! ! Purpose ! ======= ! ! DPBTF2 computes the Cholesky factorization of a real symmetric ! positive definite band matrix A. ! ! The factorization has the form ! A = U' * U , if UPLO = 'U', or ! A = L * L', if UPLO = 'L', ! where U is an upper triangular matrix, U' is the transpose of U, and ! L is lower triangular. ! ! This is the unblocked version of the algorithm, calling Level 2 BLAS. ! ! Arguments ! ========= ! ! UPLO (input) CHARACTER*1 ! Specifies whether the upper or lower triangular part of the ! symmetric matrix A is stored: ! = 'U': Upper triangular ! = 'L': Lower triangular ! ! N (input) INTEGER ! The order of the matrix A. N >= 0. ! ! KD (input) INTEGER ! The number of super-diagonals of the matrix A if UPLO = 'U', ! or the number of sub-diagonals if UPLO = 'L'. KD >= 0. ! ! AB (input/output) DOUBLE PRECISION array, dimension (LDAB,N) ! On entry, the upper or lower triangle of the symmetric band ! matrix A, stored in the first KD+1 rows of the array. The ! j-th column of A is stored in the j-th column of the array AB ! as follows: ! if UPLO = 'U', AB(kd+1+i-j,j) = A(i,j) for max(1,j-kd)<=i<=j; ! if UPLO = 'L', AB(1+i-j,j) = A(i,j) for j<=i<=min(n,j+kd). ! ! On exit, if INFO = 0, the triangular factor U or L from the ! Cholesky factorization A = U'*U or A = L*L' of the band ! matrix A, in the same storage format as A. ! ! LDAB (input) INTEGER ! The leading dimension of the array AB. LDAB >= KD+1. ! ! INFO (output) INTEGER ! = 0: successful exit ! < 0: if INFO = -k, the k-th argument had an illegal value ! > 0: if INFO = k, the leading minor of order k is not ! positive definite, and the factorization could not be ! completed. ! ! Further Details ! =============== ! ! The band storage scheme is illustrated by the following example, when ! N = 6, KD = 2, and UPLO = 'U': ! ! On entry: On exit: ! ! * * a13 a24 a35 a46 * * u13 u24 u35 u46 ! * a12 a23 a34 a45 a56 * u12 u23 u34 u45 u56 ! a11 a22 a33 a44 a55 a66 u11 u22 u33 u44 u55 u66 ! ! Similarly, if UPLO = 'L' the format of A is as follows: ! ! On entry: On exit: ! ! a11 a22 a33 a44 a55 a66 l11 l22 l33 l44 l55 l66 ! a21 a32 a43 a54 a65 * l21 l32 l43 l54 l65 * ! a31 a42 a53 a64 * * l31 l42 l53 l64 * * ! ! Array elements marked * are not used by the routine. ! ! ===================================================================== ! ! .. Parameters .. DOUBLE PRECISION ONE, ZERO PARAMETER ( ONE = 1.0D+0, ZERO = 0.5D-16 ) ! .. ! .. Local Scalars .. LOGICAL UPPER INTEGER J, KLD, KN DOUBLE PRECISION AJJ ! .. ! .. External Functions .. LOGICAL LSAME EXTERNAL LSAME ! .. ! .. External Subroutines .. EXTERNAL dscal, DSYR, XERBLA ! .. ! .. Intrinsic Functions .. INTRINSIC MAX, MIN, SQRT ! .. ! .. Executable Statements .. ! ! Test the input parameters. ! INFO = 0 UPPER = LSAME( UPLO, 'U' ) IF( .NOT.UPPER .AND. .NOT.LSAME( UPLO, 'L' ) ) THEN INFO = -1 ELSE IF( N.LT.0 ) THEN INFO = -2 ELSE IF( KD.LT.0 ) THEN INFO = -3 ELSE IF( LDAB.LT.KD+1 ) THEN INFO = -5 END IF IF( INFO.NE.0 ) THEN CALL XERBLA( 'DPBTF2', -INFO ) RETURN END IF ! ! Quick return if possible ! IF( N.EQ.0 )& RETURN ! KLD = MAX( 1, LDAB-1 ) ! IF( UPPER ) THEN ! ! Compute the Cholesky factorization A = U'*U. ! DO 10 J = 1, N ! ! Compute U(J,J) and test for non-positive-definiteness. ! AJJ = AB( KD+1, J ) IF( AJJ.LE.ZERO )& GO TO 30 AJJ = SQRT( AJJ ) AB( KD+1, J ) = AJJ ! ! Compute elements J+1:J+KN of row J and update the ! trailing submatrix within the band. ! KN = MIN( KD, N-J ) IF( KN.GT.0 ) THEN CALL dscal( KN, ONE / AJJ, AB( KD, J+1 ), KLD ) CALL DSYR( 'Upper', KN, -ONE, AB( KD, J+1 ), KLD,& AB( KD+1, J+1 ), KLD ) END IF 10 CONTINUE ELSE ! ! Compute the Cholesky factorization A = L*L'. ! DO 20 J = 1, N ! ! Compute L(J,J) and test for non-positive-definiteness. ! AJJ = AB( 1, J ) IF( AJJ.LE.ZERO )& GO TO 30 AJJ = SQRT( AJJ ) AB( 1, J ) = AJJ ! ! Compute elements J+1:J+KN of column J and update the ! trailing submatrix within the band. ! KN = MIN( KD, N-J ) IF( KN.GT.0 ) THEN CALL dscal( KN, ONE / AJJ, AB( 2, J ), 1 ) CALL DSYR( 'Lower', KN, -ONE, AB( 2, J ), 1,& AB( 1, J+1 ), KLD ) END IF 20 CONTINUE END IF RETURN ! 30 CONTINUE INFO = J RETURN ! ! End of DPBTF2 ! END SUBROUTINE DPOTF2( UPLO, N, A, LDA, INFO ) ! ! -- LAPACK routine (version 1.0) -- ! Univ. of Tennessee, Univ. of California Berkeley, NAG Ltd., ! Courant Institute, Argonne National Lab, and Rice University ! February 29, 1992 ! ! .. Scalar Arguments .. CHARACTER UPLO INTEGER INFO, LDA, N ! .. ! .. Array Arguments .. DOUBLE PRECISION A( LDA, * ) ! .. ! ! Purpose ! ======= ! ! DPOTF2 computes the Cholesky factorization of a real symmetric ! positive definite matrix A. ! ! The factorization has the form ! A = U' * U , if UPLO = 'U', or ! A = L * L', if UPLO = 'L', ! where U is an upper triangular matrix and L is lower triangular. ! ! This is the unblocked version of the algorithm, calling Level 2 BLAS. ! ! Arguments ! ========= ! ! UPLO (input) CHARACTER*1 ! Specifies whether the upper or lower triangular part of the ! symmetric matrix A is stored. ! = 'U': Upper triangular ! = 'L': Lower triangular ! ! N (input) INTEGER ! The order of the matrix A. N >= 0. ! ! A (input/output) DOUBLE PRECISION array, dimension (LDA,N) ! On entry, the symmetric matrix A. If UPLO = 'U', the leading ! n by n upper triangular part of A contains the upper ! triangular part of the matrix A, and the strictly lower ! triangular part of A is not referenced. If UPLO = 'L', the ! leading n by n lower triangular part of A contains the lower ! triangular part of the matrix A, and the strictly upper ! triangular part of A is not referenced. ! ! On exit, if INFO = 0, the factor U or L from the Cholesky ! factorization A = U'*U or A = L*L'. ! ! LDA (input) INTEGER ! The leading dimension of the array A. LDA >= max(1,N). ! ! INFO (output) INTEGER ! = 0: successful exit ! < 0: if INFO = -k, the k-th argument had an illegal value ! > 0: if INFO = k, the leading minor of order k is not ! positive definite, and the factorization could not be ! completed. ! ! ===================================================================== ! ! .. Parameters .. DOUBLE PRECISION ONE, ZERO PARAMETER ( ONE = 1.0D+0, ZERO = 0.5D-16 ) ! .. ! .. Local Scalars .. LOGICAL UPPER INTEGER J DOUBLE PRECISION AJJ ! .. ! .. External Functions .. LOGICAL LSAME DOUBLE PRECISION ddot EXTERNAL LSAME, ddot ! .. ! .. External Subroutines .. EXTERNAL DGEMV, dscal, XERBLA ! .. ! .. Intrinsic Functions .. INTRINSIC MAX, SQRT ! .. ! .. Executable Statements .. ! ! Test the input parameters. ! INFO = 0 UPPER = LSAME( UPLO, 'U' ) IF( .NOT.UPPER .AND. .NOT.LSAME( UPLO, 'L' ) ) THEN INFO = -1 ELSE IF( N.LT.0 ) THEN INFO = -2 ELSE IF( LDA.LT.MAX( 1, N ) ) THEN INFO = -4 END IF IF( INFO.NE.0 ) THEN CALL XERBLA( 'DPOTF2', -INFO ) RETURN END IF ! ! Quick return if possible ! IF( N.EQ.0 )& RETURN ! IF( UPPER ) THEN ! ! Compute the Cholesky factorization A = U'*U. ! DO 10 J = 1, N ! ! Compute U(J,J) and test for non-positive-definiteness. ! AJJ = A( J, J ) - ddot( J-1, A( 1, J ), 1, A( 1, J ), 1 ) IF( AJJ.LE.ZERO ) THEN A( J, J ) = AJJ GO TO 30 END IF AJJ = SQRT( AJJ ) A( J, J ) = AJJ ! ! Compute elements J+1:N of row J. ! IF( J.LT.N ) THEN CALL DGEMV( 'Transpose', J-1, N-J, -ONE, A( 1, J+1 ),& LDA, A( 1, J ), 1, ONE, A( J, J+1 ), LDA ) CALL dscal( N-J, ONE / AJJ, A( J, J+1 ), LDA ) END IF 10 CONTINUE ELSE ! ! Compute the Cholesky factorization A = L*L'. ! DO 20 J = 1, N ! ! Compute L(J,J) and test for non-positive-definiteness. ! AJJ = A( J, J ) - ddot( J-1, A( J, 1 ), LDA, A( J, 1 ),& LDA ) IF( AJJ.LE.ZERO ) THEN A( J, J ) = AJJ GO TO 30 END IF AJJ = SQRT( AJJ ) A( J, J ) = AJJ ! ! Compute elements J+1:N of column J. ! IF( J.LT.N ) THEN CALL DGEMV( 'No transpose', N-J, J-1, -ONE, A( J+1, 1 ),& LDA, A( J, 1 ), LDA, ONE, A( J+1, J ), 1 ) CALL dscal( N-J, ONE / AJJ, A( J+1, J ), 1 ) END IF 20 CONTINUE END IF GO TO 40 ! 30 CONTINUE INFO = J ! 40 CONTINUE RETURN ! ! End of DPOTF2 ! END SUBROUTINE DTRSM ( SIDE, UPLO, TRANSA, DIAG, M, N, ALPHA, A, LDA,& B, LDB ) ! .. Scalar Arguments .. CHARACTER*1 SIDE, UPLO, TRANSA, DIAG INTEGER M, N, LDA, LDB DOUBLE PRECISION ALPHA ! .. Array Arguments .. DOUBLE PRECISION A( LDA, * ), B( LDB, * ) ! .. ! ! Purpose ! ======= ! ! DTRSM solves one of the matrix equations ! ! op( A )*X = alpha*B, or X*op( A ) = alpha*B, ! ! where alpha is a scalar, X and B are m by n matrices, A is a unit, or ! non-unit, upper or lower triangular matrix and op( A ) is one of ! ! op( A ) = A or op( A ) = A'. ! ! The matrix X is overwritten on B. ! ! Parameters ! ========== ! ! SIDE - CHARACTER*1. ! On entry, SIDE specifies whether op( A ) appears on the left ! or right of X as follows: ! ! SIDE = 'L' or 'l' op( A )*X = alpha*B. ! ! SIDE = 'R' or 'r' X*op( A ) = alpha*B. ! ! Unchanged on exit. ! ! UPLO - CHARACTER*1. ! On entry, UPLO specifies whether the matrix A is an upper or ! lower triangular matrix as follows: ! ! UPLO = 'U' or 'u' A is an upper triangular matrix. ! ! UPLO = 'L' or 'l' A is a lower triangular matrix. ! ! Unchanged on exit. ! ! TRANSA - CHARACTER*1. ! On entry, TRANSA specifies the form of op( A ) to be used in ! the matrix multiplication as follows: ! ! TRANSA = 'N' or 'n' op( A ) = A. ! ! TRANSA = 'T' or 't' op( A ) = A'. ! ! TRANSA = 'C' or 'c' op( A ) = A'. ! ! Unchanged on exit. ! ! DIAG - CHARACTER*1. ! On entry, DIAG specifies whether or not A is unit triangular ! as follows: ! ! DIAG = 'U' or 'u' A is assumed to be unit triangular. ! ! DIAG = 'N' or 'n' A is not assumed to be unit ! triangular. ! ! Unchanged on exit. ! ! M - INTEGER. ! On entry, M specifies the number of rows of B. M must be at ! least zero. ! Unchanged on exit. ! ! N - INTEGER. ! On entry, N specifies the number of columns of B. N must be ! at least zero. ! Unchanged on exit. ! ! ALPHA - DOUBLE PRECISION. ! On entry, ALPHA specifies the scalar alpha. When alpha is ! zero then A is not referenced and B need not be set before ! entry. ! Unchanged on exit. ! ! A - DOUBLE PRECISION array of DIMENSION ( LDA, k ), where k is m ! when SIDE = 'L' or 'l' and is n when SIDE = 'R' or 'r'. ! Before entry with UPLO = 'U' or 'u', the leading k by k ! upper triangular part of the array A must contain the upper ! triangular matrix and the strictly lower triangular part of ! A is not referenced. ! Before entry with UPLO = 'L' or 'l', the leading k by k ! lower triangular part of the array A must contain the lower ! triangular matrix and the strictly upper triangular part of ! A is not referenced. ! Note that when DIAG = 'U' or 'u', the diagonal elements of ! A are not referenced either, but are assumed to be unity. ! Unchanged on exit. ! ! LDA - INTEGER. ! On entry, LDA specifies the first dimension of A as declared ! in the calling (sub) program. When SIDE = 'L' or 'l' then ! LDA must be at least max( 1, m ), when SIDE = 'R' or 'r' ! then LDA must be at least max( 1, n ). ! Unchanged on exit. ! ! B - DOUBLE PRECISION array of DIMENSION ( LDB, n ). ! Before entry, the leading m by n part of the array B must ! contain the right-hand side matrix B, and on exit is ! overwritten by the solution matrix X. ! ! LDB - INTEGER. ! On entry, LDB specifies the first dimension of B as declared ! in the calling (sub) program. LDB must be at least ! max( 1, m ). ! Unchanged on exit. ! ! ! Level 3 Blas routine. ! ! ! -- Written on 8-February-1989. ! Jack Dongarra, Argonne National Laboratory. ! Iain Duff, AERE Harwell. ! Jeremy Du Croz, Numerical Algorithms Group Ltd. ! Sven Hammarling, Numerical Algorithms Group Ltd. ! ! ! .. External Functions .. LOGICAL LSAME EXTERNAL LSAME ! .. External Subroutines .. EXTERNAL XERBLA ! .. Intrinsic Functions .. INTRINSIC MAX ! .. Local Scalars .. LOGICAL LSIDE, NOUNIT, UPPER INTEGER I, INFO, J, K, NROWA DOUBLE PRECISION TEMP ! .. Parameters .. DOUBLE PRECISION ONE , ZERO PARAMETER ( ONE = 1.0D+0, ZERO = 0.0D+0 ) ! .. ! .. Executable Statements .. ! ! Test the input parameters. ! LSIDE = LSAME( SIDE , 'L' ) IF( LSIDE )THEN NROWA = M ELSE NROWA = N END IF NOUNIT = LSAME( DIAG , 'N' ) UPPER = LSAME( UPLO , 'U' ) ! INFO = 0 IF( ( .NOT.LSIDE ).AND.& ( .NOT.LSAME( SIDE , 'R' ) ) )THEN INFO = 1 ELSE IF( ( .NOT.UPPER ).AND.& ( .NOT.LSAME( UPLO , 'L' ) ) )THEN INFO = 2 ELSE IF( ( .NOT.LSAME( TRANSA, 'N' ) ).AND.& ( .NOT.LSAME( TRANSA, 'T' ) ).AND.& ( .NOT.LSAME( TRANSA, 'C' ) ) )THEN INFO = 3 ELSE IF( ( .NOT.LSAME( DIAG , 'U' ) ).AND.& ( .NOT.LSAME( DIAG , 'N' ) ) )THEN INFO = 4 ELSE IF( M .LT.0 )THEN INFO = 5 ELSE IF( N .LT.0 )THEN INFO = 6 ELSE IF( LDA.LT.MAX( 1, NROWA ) )THEN INFO = 9 ELSE IF( LDB.LT.MAX( 1, M ) )THEN INFO = 11 END IF IF( INFO.NE.0 )THEN CALL XERBLA( 'DTRSM ', INFO ) RETURN END IF ! ! Quick return if possible. ! IF( N.EQ.0 )& RETURN ! ! And when alpha.eq.zero. ! IF( ALPHA.EQ.ZERO )THEN DO 20, J = 1, N DO 10, I = 1, M B( I, J ) = ZERO 10 CONTINUE 20 CONTINUE RETURN END IF ! ! Start the operations. ! IF( LSIDE )THEN IF( LSAME( TRANSA, 'N' ) )THEN ! ! Form B := alpha*inv( A )*B. ! IF( UPPER )THEN DO 60, J = 1, N IF( ALPHA.NE.ONE )THEN DO 30, I = 1, M B( I, J ) = ALPHA*B( I, J ) 30 CONTINUE END IF DO 50, K = M, 1, -1 IF( B( K, J ).NE.ZERO )THEN IF( NOUNIT )& B( K, J ) = B( K, J )/A( K, K ) DO 40, I = 1, K - 1 B( I, J ) = B( I, J ) - B( K, J )*A( I, K ) 40 CONTINUE END IF 50 CONTINUE 60 CONTINUE ELSE DO 100, J = 1, N IF( ALPHA.NE.ONE )THEN DO 70, I = 1, M B( I, J ) = ALPHA*B( I, J ) 70 CONTINUE END IF DO 90 K = 1, M IF( B( K, J ).NE.ZERO )THEN IF( NOUNIT )& B( K, J ) = B( K, J )/A( K, K ) DO 80, I = K + 1, M B( I, J ) = B( I, J ) - B( K, J )*A( I, K ) 80 CONTINUE END IF 90 CONTINUE 100 CONTINUE END IF ELSE ! ! Form B := alpha*inv( A' )*B. ! IF( UPPER )THEN DO 130, J = 1, N DO 120, I = 1, M TEMP = ALPHA*B( I, J ) DO 110, K = 1, I - 1 TEMP = TEMP - A( K, I )*B( K, J ) 110 CONTINUE IF( NOUNIT )& TEMP = TEMP/A( I, I ) B( I, J ) = TEMP 120 CONTINUE 130 CONTINUE ELSE DO 160, J = 1, N DO 150, I = M, 1, -1 TEMP = ALPHA*B( I, J ) DO 140, K = I + 1, M TEMP = TEMP - A( K, I )*B( K, J ) 140 CONTINUE IF( NOUNIT )& TEMP = TEMP/A( I, I ) B( I, J ) = TEMP 150 CONTINUE 160 CONTINUE END IF END IF ELSE IF( LSAME( TRANSA, 'N' ) )THEN ! ! Form B := alpha*B*inv( A ). ! IF( UPPER )THEN DO 210, J = 1, N IF( ALPHA.NE.ONE )THEN DO 170, I = 1, M B( I, J ) = ALPHA*B( I, J ) 170 CONTINUE END IF DO 190, K = 1, J - 1 IF( A( K, J ).NE.ZERO )THEN DO 180, I = 1, M B( I, J ) = B( I, J ) - A( K, J )*B( I, K ) 180 CONTINUE END IF 190 CONTINUE IF( NOUNIT )THEN TEMP = ONE/A( J, J ) DO 200, I = 1, M B( I, J ) = TEMP*B( I, J ) 200 CONTINUE END IF 210 CONTINUE ELSE DO 260, J = N, 1, -1 IF( ALPHA.NE.ONE )THEN DO 220, I = 1, M B( I, J ) = ALPHA*B( I, J ) 220 CONTINUE END IF DO 240, K = J + 1, N IF( A( K, J ).NE.ZERO )THEN DO 230, I = 1, M B( I, J ) = B( I, J ) - A( K, J )*B( I, K ) 230 CONTINUE END IF 240 CONTINUE IF( NOUNIT )THEN TEMP = ONE/A( J, J ) DO 250, I = 1, M B( I, J ) = TEMP*B( I, J ) 250 CONTINUE END IF 260 CONTINUE END IF ELSE ! ! Form B := alpha*B*inv( A' ). ! IF( UPPER )THEN DO 310, K = N, 1, -1 IF( NOUNIT )THEN TEMP = ONE/A( K, K ) DO 270, I = 1, M B( I, K ) = TEMP*B( I, K ) 270 CONTINUE END IF DO 290, J = 1, K - 1 IF( A( J, K ).NE.ZERO )THEN TEMP = A( J, K ) DO 280, I = 1, M B( I, J ) = B( I, J ) - TEMP*B( I, K ) 280 CONTINUE END IF 290 CONTINUE IF( ALPHA.NE.ONE )THEN DO 300, I = 1, M B( I, K ) = ALPHA*B( I, K ) 300 CONTINUE END IF 310 CONTINUE ELSE DO 360, K = 1, N IF( NOUNIT )THEN TEMP = ONE/A( K, K ) DO 320, I = 1, M B( I, K ) = TEMP*B( I, K ) 320 CONTINUE END IF DO 340, J = K + 1, N IF( A( J, K ).NE.ZERO )THEN TEMP = A( J, K ) DO 330, I = 1, M B( I, J ) = B( I, J ) - TEMP*B( I, K ) 330 CONTINUE END IF 340 CONTINUE IF( ALPHA.NE.ONE )THEN DO 350, I = 1, M B( I, K ) = ALPHA*B( I, K ) 350 CONTINUE END IF 360 CONTINUE END IF END IF END IF ! RETURN ! ! End of DTRSM . ! END SUBROUTINE DSYRK ( UPLO, TRANS, N, K, ALPHA, A, LDA,& BETA, C, LDC ) ! .. Scalar Arguments .. CHARACTER*1 UPLO, TRANS INTEGER N, K, LDA, LDC DOUBLE PRECISION ALPHA, BETA ! .. Array Arguments .. DOUBLE PRECISION A( LDA, * ), C( LDC, * ) ! .. ! ! Purpose ! ======= ! ! DSYRK performs one of the symmetric rank k operations ! ! C := alpha*A*A' + beta*C, ! ! or ! ! C := alpha*A'*A + beta*C, ! ! where alpha and beta are scalars, C is an n by n symmetric matrix ! and A is an n by k matrix in the first case and a k by n matrix ! in the second case. ! ! Parameters ! ========== ! ! UPLO - CHARACTER*1. ! On entry, UPLO specifies whether the upper or lower ! triangular part of the array C is to be referenced as ! follows: ! ! UPLO = 'U' or 'u' Only the upper triangular part of C ! is to be referenced. ! ! UPLO = 'L' or 'l' Only the lower triangular part of C ! is to be referenced. ! ! Unchanged on exit. ! ! TRANS - CHARACTER*1. ! On entry, TRANS specifies the operation to be performed as ! follows: ! ! TRANS = 'N' or 'n' C := alpha*A*A' + beta*C. ! ! TRANS = 'T' or 't' C := alpha*A'*A + beta*C. ! ! TRANS = 'C' or 'c' C := alpha*A'*A + beta*C. ! ! Unchanged on exit. ! ! N - INTEGER. ! On entry, N specifies the order of the matrix C. N must be ! at least zero. ! Unchanged on exit. ! ! K - INTEGER. ! On entry with TRANS = 'N' or 'n', K specifies the number ! of columns of the matrix A, and on entry with ! TRANS = 'T' or 't' or 'C' or 'c', K specifies the number ! of rows of the matrix A. K must be at least zero. ! Unchanged on exit. ! ! ALPHA - DOUBLE PRECISION. ! On entry, ALPHA specifies the scalar alpha. ! Unchanged on exit. ! ! A - DOUBLE PRECISION array of DIMENSION ( LDA, ka ), where ka is ! k when TRANS = 'N' or 'n', and is n otherwise. ! Before entry with TRANS = 'N' or 'n', the leading n by k ! part of the array A must contain the matrix A, otherwise ! the leading k by n part of the array A must contain the ! matrix A. ! Unchanged on exit. ! ! LDA - INTEGER. ! On entry, LDA specifies the first dimension of A as declared ! in the calling (sub) program. When TRANS = 'N' or 'n' ! then LDA must be at least max( 1, n ), otherwise LDA must ! be at least max( 1, k ). ! Unchanged on exit. ! ! BETA - DOUBLE PRECISION. ! On entry, BETA specifies the scalar beta. ! Unchanged on exit. ! ! C - DOUBLE PRECISION array of DIMENSION ( LDC, n ). ! Before entry with UPLO = 'U' or 'u', the leading n by n ! upper triangular part of the array C must contain the upper ! triangular part of the symmetric matrix and the strictly ! lower triangular part of C is not referenced. On exit, the ! upper triangular part of the array C is overwritten by the ! upper triangular part of the updated matrix. ! Before entry with UPLO = 'L' or 'l', the leading n by n ! lower triangular part of the array C must contain the lower ! triangular part of the symmetric matrix and the strictly ! upper triangular part of C is not referenced. On exit, the ! lower triangular part of the array C is overwritten by the ! lower triangular part of the updated matrix. ! ! LDC - INTEGER. ! On entry, LDC specifies the first dimension of C as declared ! in the calling (sub) program. LDC must be at least ! max( 1, n ). ! Unchanged on exit. ! ! ! Level 3 Blas routine. ! ! -- Written on 8-February-1989. ! Jack Dongarra, Argonne National Laboratory. ! Iain Duff, AERE Harwell. ! Jeremy Du Croz, Numerical Algorithms Group Ltd. ! Sven Hammarling, Numerical Algorithms Group Ltd. ! ! ! .. External Functions .. LOGICAL LSAME EXTERNAL LSAME ! .. External Subroutines .. EXTERNAL XERBLA ! .. Intrinsic Functions .. INTRINSIC MAX ! .. Local Scalars .. LOGICAL UPPER INTEGER I, INFO, J, L, NROWA DOUBLE PRECISION TEMP ! .. Parameters .. DOUBLE PRECISION ONE , ZERO PARAMETER ( ONE = 1.0D+0, ZERO = 0.0D+0 ) ! .. ! .. Executable Statements .. ! ! Test the input parameters. ! IF( LSAME( TRANS, 'N' ) )THEN NROWA = N ELSE NROWA = K END IF UPPER = LSAME( UPLO, 'U' ) ! INFO = 0 IF( ( .NOT.UPPER ).AND.& ( .NOT.LSAME( UPLO , 'L' ) ) )THEN INFO = 1 ELSE IF( ( .NOT.LSAME( TRANS, 'N' ) ).AND.& ( .NOT.LSAME( TRANS, 'T' ) ).AND.& ( .NOT.LSAME( TRANS, 'C' ) ) )THEN INFO = 2 ELSE IF( N .LT.0 )THEN INFO = 3 ELSE IF( K .LT.0 )THEN INFO = 4 ELSE IF( LDA.LT.MAX( 1, NROWA ) )THEN INFO = 7 ELSE IF( LDC.LT.MAX( 1, N ) )THEN INFO = 10 END IF IF( INFO.NE.0 )THEN CALL XERBLA( 'DSYRK ', INFO ) RETURN END IF ! ! Quick return if possible. ! IF( ( N.EQ.0 ).OR.& ( ( ( ALPHA.EQ.ZERO ).OR.( K.EQ.0 ) ).AND.( BETA.EQ.ONE ) ) )& RETURN ! ! And when alpha.eq.zero. ! IF( ALPHA.EQ.ZERO )THEN IF( UPPER )THEN IF( BETA.EQ.ZERO )THEN DO 20, J = 1, N DO 10, I = 1, J C( I, J ) = ZERO 10 CONTINUE 20 CONTINUE ELSE DO 40, J = 1, N DO 30, I = 1, J C( I, J ) = BETA*C( I, J ) 30 CONTINUE 40 CONTINUE END IF ELSE IF( BETA.EQ.ZERO )THEN DO 60, J = 1, N DO 50, I = J, N C( I, J ) = ZERO 50 CONTINUE 60 CONTINUE ELSE DO 80, J = 1, N DO 70, I = J, N C( I, J ) = BETA*C( I, J ) 70 CONTINUE 80 CONTINUE END IF END IF RETURN END IF ! ! Start the operations. ! IF( LSAME( TRANS, 'N' ) )THEN ! ! Form C := alpha*A*A' + beta*C. ! IF( UPPER )THEN DO 130, J = 1, N IF( BETA.EQ.ZERO )THEN DO 90, I = 1, J C( I, J ) = ZERO 90 CONTINUE ELSE IF( BETA.NE.ONE )THEN DO 100, I = 1, J C( I, J ) = BETA*C( I, J ) 100 CONTINUE END IF DO 120, L = 1, K IF( A( J, L ).NE.ZERO )THEN TEMP = ALPHA*A( J, L ) DO 110, I = 1, J C( I, J ) = C( I, J ) + TEMP*A( I, L ) 110 CONTINUE END IF 120 CONTINUE 130 CONTINUE ELSE DO 180, J = 1, N IF( BETA.EQ.ZERO )THEN DO 140, I = J, N C( I, J ) = ZERO 140 CONTINUE ELSE IF( BETA.NE.ONE )THEN DO 150, I = J, N C( I, J ) = BETA*C( I, J ) 150 CONTINUE END IF DO 170, L = 1, K IF( A( J, L ).NE.ZERO )THEN TEMP = ALPHA*A( J, L ) DO 160, I = J, N C( I, J ) = C( I, J ) + TEMP*A( I, L ) 160 CONTINUE END IF 170 CONTINUE 180 CONTINUE END IF ELSE ! ! Form C := alpha*A'*A + beta*C. ! IF( UPPER )THEN DO 210, J = 1, N DO 200, I = 1, J TEMP = ZERO DO 190, L = 1, K TEMP = TEMP + A( L, I )*A( L, J ) 190 CONTINUE IF( BETA.EQ.ZERO )THEN C( I, J ) = ALPHA*TEMP ELSE C( I, J ) = ALPHA*TEMP + BETA*C( I, J ) END IF 200 CONTINUE 210 CONTINUE ELSE DO 240, J = 1, N DO 230, I = J, N TEMP = ZERO DO 220, L = 1, K TEMP = TEMP + A( L, I )*A( L, J ) 220 CONTINUE IF( BETA.EQ.ZERO )THEN C( I, J ) = ALPHA*TEMP ELSE C( I, J ) = ALPHA*TEMP + BETA*C( I, J ) END IF 230 CONTINUE 240 CONTINUE END IF END IF ! RETURN ! ! End of DSYRK . ! END SUBROUTINE DGEMM ( TRANSA, TRANSB, M, N, K, ALPHA, A, LDA, B, LDB,& BETA, C, LDC ) ! .. Scalar Arguments .. CHARACTER*1 TRANSA, TRANSB INTEGER M, N, K, LDA, LDB, LDC DOUBLE PRECISION ALPHA, BETA ! .. Array Arguments .. DOUBLE PRECISION A( LDA, * ), B( LDB, * ), C( LDC, * ) ! .. ! ! Purpose ! ======= ! ! DGEMM performs one of the matrix-matrix operations ! ! C := alpha*op( A )*op( B ) + beta*C, ! ! where op( X ) is one of ! ! op( X ) = X or op( X ) = X', ! ! alpha and beta are scalars, and A, B and C are matrices, with op( A ) ! an m by k matrix, op( B ) a k by n matrix and C an m by n matrix. ! ! Parameters ! ========== ! ! TRANSA - CHARACTER*1. ! On entry, TRANSA specifies the form of op( A ) to be used in ! the matrix multiplication as follows: ! ! TRANSA = 'N' or 'n', op( A ) = A. ! ! TRANSA = 'T' or 't', op( A ) = A'. ! ! TRANSA = 'C' or 'c', op( A ) = A'. ! ! Unchanged on exit. ! ! TRANSB - CHARACTER*1. ! On entry, TRANSB specifies the form of op( B ) to be used in ! the matrix multiplication as follows: ! ! TRANSB = 'N' or 'n', op( B ) = B. ! ! TRANSB = 'T' or 't', op( B ) = B'. ! ! TRANSB = 'C' or 'c', op( B ) = B'. ! ! Unchanged on exit. ! ! M - INTEGER. ! On entry, M specifies the number of rows of the matrix ! op( A ) and of the matrix C. M must be at least zero. ! Unchanged on exit. ! ! N - INTEGER. ! On entry, N specifies the number of columns of the matrix ! op( B ) and the number of columns of the matrix C. N must be ! at least zero. ! Unchanged on exit. ! ! K - INTEGER. ! On entry, K specifies the number of columns of the matrix ! op( A ) and the number of rows of the matrix op( B ). K must ! be at least zero. ! Unchanged on exit. ! ! ALPHA - DOUBLE PRECISION. ! On entry, ALPHA specifies the scalar alpha. ! Unchanged on exit. ! ! A - DOUBLE PRECISION array of DIMENSION ( LDA, ka ), where ka is ! k when TRANSA = 'N' or 'n', and is m otherwise. ! Before entry with TRANSA = 'N' or 'n', the leading m by k ! part of the array A must contain the matrix A, otherwise ! the leading k by m part of the array A must contain the ! matrix A. ! Unchanged on exit. ! ! LDA - INTEGER. ! On entry, LDA specifies the first dimension of A as declared ! in the calling (sub) program. When TRANSA = 'N' or 'n' then ! LDA must be at least max( 1, m ), otherwise LDA must be at ! least max( 1, k ). ! Unchanged on exit. ! ! B - DOUBLE PRECISION array of DIMENSION ( LDB, kb ), where kb is ! n when TRANSB = 'N' or 'n', and is k otherwise. ! Before entry with TRANSB = 'N' or 'n', the leading k by n ! part of the array B must contain the matrix B, otherwise ! the leading n by k part of the array B must contain the ! matrix B. ! Unchanged on exit. ! ! LDB - INTEGER. ! On entry, LDB specifies the first dimension of B as declared ! in the calling (sub) program. When TRANSB = 'N' or 'n' then ! LDB must be at least max( 1, k ), otherwise LDB must be at ! least max( 1, n ). ! Unchanged on exit. ! ! BETA - DOUBLE PRECISION. ! On entry, BETA specifies the scalar beta. When BETA is ! supplied as zero then C need not be set on input. ! Unchanged on exit. ! ! C - DOUBLE PRECISION array of DIMENSION ( LDC, n ). ! Before entry, the leading m by n part of the array C must ! contain the matrix C, except when beta is zero, in which ! case C need not be set on entry. ! On exit, the array C is overwritten by the m by n matrix ! ( alpha*op( A )*op( B ) + beta*C ). ! ! LDC - INTEGER. ! On entry, LDC specifies the first dimension of C as declared ! in the calling (sub) program. LDC must be at least ! max( 1, m ). ! Unchanged on exit. ! ! ! Level 3 Blas routine. ! ! -- Written on 8-February-1989. ! Jack Dongarra, Argonne National Laboratory. ! Iain Duff, AERE Harwell. ! Jeremy Du Croz, Numerical Algorithms Group Ltd. ! Sven Hammarling, Numerical Algorithms Group Ltd. ! ! ! .. External Functions .. LOGICAL LSAME EXTERNAL LSAME ! .. External Subroutines .. EXTERNAL XERBLA ! .. Intrinsic Functions .. INTRINSIC MAX ! .. Local Scalars .. LOGICAL NOTA, NOTB INTEGER I, INFO, J, L, NCOLA, NROWA, NROWB DOUBLE PRECISION TEMP ! .. Parameters .. DOUBLE PRECISION ONE , ZERO PARAMETER ( ONE = 1.0D+0, ZERO = 0.0D+0 ) ! .. ! .. Executable Statements .. ! ! Set NOTA and NOTB as true if A and B respectively are not ! transposed and set NROWA, NCOLA and NROWB as the number of rows ! and columns of A and the number of rows of B respectively. ! NOTA = LSAME( TRANSA, 'N' ) NOTB = LSAME( TRANSB, 'N' ) IF( NOTA )THEN NROWA = M NCOLA = K ELSE NROWA = K NCOLA = M END IF IF( NOTB )THEN NROWB = K ELSE NROWB = N END IF ! ! Test the input parameters. ! INFO = 0 IF( ( .NOT.NOTA ).AND.& ( .NOT.LSAME( TRANSA, 'C' ) ).AND.& ( .NOT.LSAME( TRANSA, 'T' ) ) )THEN INFO = 1 ELSE IF( ( .NOT.NOTB ).AND.& ( .NOT.LSAME( TRANSB, 'C' ) ).AND.& ( .NOT.LSAME( TRANSB, 'T' ) ) )THEN INFO = 2 ELSE IF( M .LT.0 )THEN INFO = 3 ELSE IF( N .LT.0 )THEN INFO = 4 ELSE IF( K .LT.0 )THEN INFO = 5 ELSE IF( LDA.LT.MAX( 1, NROWA ) )THEN INFO = 8 ELSE IF( LDB.LT.MAX( 1, NROWB ) )THEN INFO = 10 ELSE IF( LDC.LT.MAX( 1, M ) )THEN INFO = 13 END IF IF( INFO.NE.0 )THEN CALL XERBLA( 'DGEMM ', INFO ) RETURN END IF ! ! Quick return if possible. ! IF( ( M.EQ.0 ).OR.( N.EQ.0 ).OR.& ( ( ( ALPHA.EQ.ZERO ).OR.( K.EQ.0 ) ).AND.( BETA.EQ.ONE ) ) )& RETURN ! ! And if alpha.eq.zero. ! IF( ALPHA.EQ.ZERO )THEN IF( BETA.EQ.ZERO )THEN DO 20, J = 1, N DO 10, I = 1, M C( I, J ) = ZERO 10 CONTINUE 20 CONTINUE ELSE DO 40, J = 1, N DO 30, I = 1, M C( I, J ) = BETA*C( I, J ) 30 CONTINUE 40 CONTINUE END IF RETURN END IF ! ! Start the operations. ! IF( NOTB )THEN IF( NOTA )THEN ! ! Form C := alpha*A*B + beta*C. ! DO 90, J = 1, N IF( BETA.EQ.ZERO )THEN DO 50, I = 1, M C( I, J ) = ZERO 50 CONTINUE ELSE IF( BETA.NE.ONE )THEN DO 60, I = 1, M C( I, J ) = BETA*C( I, J ) 60 CONTINUE END IF DO 80, L = 1, K IF( B( L, J ).NE.ZERO )THEN TEMP = ALPHA*B( L, J ) DO 70, I = 1, M C( I, J ) = C( I, J ) + TEMP*A( I, L ) 70 CONTINUE END IF 80 CONTINUE 90 CONTINUE ELSE ! ! Form C := alpha*A'*B + beta*C ! DO 120, J = 1, N DO 110, I = 1, M TEMP = ZERO DO 100, L = 1, K TEMP = TEMP + A( L, I )*B( L, J ) 100 CONTINUE IF( BETA.EQ.ZERO )THEN C( I, J ) = ALPHA*TEMP ELSE C( I, J ) = ALPHA*TEMP + BETA*C( I, J ) END IF 110 CONTINUE 120 CONTINUE END IF ELSE IF( NOTA )THEN ! ! Form C := alpha*A*B' + beta*C ! DO 170, J = 1, N IF( BETA.EQ.ZERO )THEN DO 130, I = 1, M C( I, J ) = ZERO 130 CONTINUE ELSE IF( BETA.NE.ONE )THEN DO 140, I = 1, M C( I, J ) = BETA*C( I, J ) 140 CONTINUE END IF DO 160, L = 1, K IF( B( J, L ).NE.ZERO )THEN TEMP = ALPHA*B( J, L ) DO 150, I = 1, M C( I, J ) = C( I, J ) + TEMP*A( I, L ) 150 CONTINUE END IF 160 CONTINUE 170 CONTINUE ELSE ! ! Form C := alpha*A'*B' + beta*C ! DO 200, J = 1, N DO 190, I = 1, M TEMP = ZERO DO 180, L = 1, K TEMP = TEMP + A( L, I )*B( J, L ) 180 CONTINUE IF( BETA.EQ.ZERO )THEN C( I, J ) = ALPHA*TEMP ELSE C( I, J ) = ALPHA*TEMP + BETA*C( I, J ) END IF 190 CONTINUE 200 CONTINUE END IF END IF ! RETURN ! ! End of DGEMM . ! END SUBROUTINE DTBSV ( UPLO, TRANS, DIAG, N, K, A, LDA, X, INCX ) ! .. Scalar Arguments .. INTEGER INCX, K, LDA, N CHARACTER*1 DIAG, TRANS, UPLO ! .. Array Arguments .. DOUBLE PRECISION A( LDA, * ), X( * ) ! .. ! ! Purpose ! ======= ! ! DTBSV solves one of the systems of equations ! ! A*x = b, or A'*x = b, ! ! where b and x are n element vectors and A is an n by n unit, or ! non-unit, upper or lower triangular band matrix, with ( k + 1 ) ! diagonals. ! ! No test for singularity or near-singularity is included in this ! routine. Such tests must be performed before calling this routine. ! ! Parameters ! ========== ! ! UPLO - CHARACTER*1. ! On entry, UPLO specifies whether the matrix is an upper or ! lower triangular matrix as follows: ! ! UPLO = 'U' or 'u' A is an upper triangular matrix. ! ! UPLO = 'L' or 'l' A is a lower triangular matrix. ! ! Unchanged on exit. ! ! TRANS - CHARACTER*1. ! On entry, TRANS specifies the equations to be solved as ! follows: ! ! TRANS = 'N' or 'n' A*x = b. ! ! TRANS = 'T' or 't' A'*x = b. ! ! TRANS = 'C' or 'c' A'*x = b. ! ! Unchanged on exit. ! ! DIAG - CHARACTER*1. ! On entry, DIAG specifies whether or not A is unit ! triangular as follows: ! ! DIAG = 'U' or 'u' A is assumed to be unit triangular. ! ! DIAG = 'N' or 'n' A is not assumed to be unit ! triangular. ! ! Unchanged on exit. ! ! N - INTEGER. ! On entry, N specifies the order of the matrix A. ! N must be at least zero. ! Unchanged on exit. ! ! K - INTEGER. ! On entry with UPLO = 'U' or 'u', K specifies the number of ! super-diagonals of the matrix A. ! On entry with UPLO = 'L' or 'l', K specifies the number of ! sub-diagonals of the matrix A. ! K must satisfy 0 .le. K. ! Unchanged on exit. ! ! A - DOUBLE PRECISION array of DIMENSION ( LDA, n ). ! Before entry with UPLO = 'U' or 'u', the leading ( k + 1 ) ! by n part of the array A must contain the upper triangular ! band part of the matrix of coefficients, supplied column by ! column, with the leading diagonal of the matrix in row ! ( k + 1 ) of the array, the first super-diagonal starting at ! position 2 in row k, and so on. The top left k by k triangle ! of the array A is not referenced. ! The following program segment will transfer an upper ! triangular band matrix from conventional full matrix storage ! to band storage: ! ! DO 20, J = 1, N ! M = K + 1 - J ! DO 10, I = MAX( 1, J - K ), J ! A( M + I, J ) = matrix( I, J ) ! 10 CONTINUE ! 20 CONTINUE ! ! Before entry with UPLO = 'L' or 'l', the leading ( k + 1 ) ! by n part of the array A must contain the lower triangular ! band part of the matrix of coefficients, supplied column by ! column, with the leading diagonal of the matrix in row 1 of ! the array, the first sub-diagonal starting at position 1 in ! row 2, and so on. The bottom right k by k triangle of the ! array A is not referenced. ! The following program segment will transfer a lower ! triangular band matrix from conventional full matrix storage ! to band storage: ! ! DO 20, J = 1, N ! M = 1 - J ! DO 10, I = J, MIN( N, J + K ) ! A( M + I, J ) = matrix( I, J ) ! 10 CONTINUE ! 20 CONTINUE ! ! Note that when DIAG = 'U' or 'u' the elements of the array A ! corresponding to the diagonal elements of the matrix are not ! referenced, but are assumed to be unity. ! Unchanged on exit. ! ! LDA - INTEGER. ! On entry, LDA specifies the first dimension of A as declared ! in the calling (sub) program. LDA must be at least ! ( k + 1 ). ! Unchanged on exit. ! ! X - DOUBLE PRECISION array of dimension at least ! ( 1 + ( n - 1 )*abs( INCX ) ). ! Before entry, the incremented array X must contain the n ! element right-hand side vector b. On exit, X is overwritten ! with the solution vector x. ! ! INCX - INTEGER. ! On entry, INCX specifies the increment for the elements of ! X. INCX must not be zero. ! Unchanged on exit. ! ! ! Level 2 Blas routine. ! ! -- Written on 22-October-1986. ! Jack Dongarra, Argonne National Lab. ! Jeremy Du Croz, Nag Central Office. ! Sven Hammarling, Nag Central Office. ! Richard Hanson, Sandia National Labs. ! ! ! .. Parameters .. DOUBLE PRECISION ZERO PARAMETER ( ZERO = 0.0D+0 ) ! .. Local Scalars .. DOUBLE PRECISION TEMP INTEGER I, INFO, IX, J, JX, KPLUS1, KX, L LOGICAL NOUNIT ! .. External Functions .. LOGICAL LSAME EXTERNAL LSAME ! .. External Subroutines .. EXTERNAL XERBLA ! .. Intrinsic Functions .. INTRINSIC MAX, MIN ! .. ! .. Executable Statements .. ! ! Test the input parameters. ! INFO = 0 IF ( .NOT.LSAME( UPLO , 'U' ).AND.& .NOT.LSAME( UPLO , 'L' ) )THEN INFO = 1 ELSE IF( .NOT.LSAME( TRANS, 'N' ).AND.& .NOT.LSAME( TRANS, 'T' ).AND.& .NOT.LSAME( TRANS, 'C' ) )THEN INFO = 2 ELSE IF( .NOT.LSAME( DIAG , 'U' ).AND.& .NOT.LSAME( DIAG , 'N' ) )THEN INFO = 3 ELSE IF( N.LT.0 )THEN INFO = 4 ELSE IF( K.LT.0 )THEN INFO = 5 ELSE IF( LDA.LT.( K + 1 ) )THEN INFO = 7 ELSE IF( INCX.EQ.0 )THEN INFO = 9 END IF IF( INFO.NE.0 )THEN CALL XERBLA( 'DTBSV ', INFO ) RETURN END IF ! ! Quick return if possible. ! IF( N.EQ.0 )& RETURN ! NOUNIT = LSAME( DIAG, 'N' ) ! ! Set up the start point in X if the increment is not unity. This ! will be ( N - 1 )*INCX too small for descending loops. ! IF( INCX.LE.0 )THEN KX = 1 - ( N - 1 )*INCX ELSE IF( INCX.NE.1 )THEN KX = 1 END IF ! ! Start the operations. In this version the elements of A are ! accessed by sequentially with one pass through A. ! IF( LSAME( TRANS, 'N' ) )THEN ! ! Form x := inv( A )*x. ! IF( LSAME( UPLO, 'U' ) )THEN KPLUS1 = K + 1 IF( INCX.EQ.1 )THEN DO 20, J = N, 1, -1 IF( X( J ).NE.ZERO )THEN L = KPLUS1 - J IF( NOUNIT )& X( J ) = X( J )/A( KPLUS1, J ) TEMP = X( J ) DO 10, I = J - 1, MAX( 1, J - K ), -1 X( I ) = X( I ) - TEMP*A( L + I, J ) 10 CONTINUE END IF 20 CONTINUE ELSE KX = KX + ( N - 1 )*INCX JX = KX DO 40, J = N, 1, -1 KX = KX - INCX IF( X( JX ).NE.ZERO )THEN IX = KX L = KPLUS1 - J IF( NOUNIT )& X( JX ) = X( JX )/A( KPLUS1, J ) TEMP = X( JX ) DO 30, I = J - 1, MAX( 1, J - K ), -1 X( IX ) = X( IX ) - TEMP*A( L + I, J ) IX = IX - INCX 30 CONTINUE END IF JX = JX - INCX 40 CONTINUE END IF ELSE IF( INCX.EQ.1 )THEN DO 60, J = 1, N IF( X( J ).NE.ZERO )THEN L = 1 - J IF( NOUNIT )& X( J ) = X( J )/A( 1, J ) TEMP = X( J ) DO 50, I = J + 1, MIN( N, J + K ) X( I ) = X( I ) - TEMP*A( L + I, J ) 50 CONTINUE END IF 60 CONTINUE ELSE JX = KX DO 80, J = 1, N KX = KX + INCX IF( X( JX ).NE.ZERO )THEN IX = KX L = 1 - J IF( NOUNIT )& X( JX ) = X( JX )/A( 1, J ) TEMP = X( JX ) DO 70, I = J + 1, MIN( N, J + K ) X( IX ) = X( IX ) - TEMP*A( L + I, J ) IX = IX + INCX 70 CONTINUE END IF JX = JX + INCX 80 CONTINUE END IF END IF ELSE ! ! Form x := inv( A')*x. ! IF( LSAME( UPLO, 'U' ) )THEN KPLUS1 = K + 1 IF( INCX.EQ.1 )THEN DO 100, J = 1, N TEMP = X( J ) L = KPLUS1 - J DO 90, I = MAX( 1, J - K ), J - 1 TEMP = TEMP - A( L + I, J )*X( I ) 90 CONTINUE IF( NOUNIT )& TEMP = TEMP/A( KPLUS1, J ) X( J ) = TEMP 100 CONTINUE ELSE JX = KX DO 120, J = 1, N TEMP = X( JX ) IX = KX L = KPLUS1 - J DO 110, I = MAX( 1, J - K ), J - 1 TEMP = TEMP - A( L + I, J )*X( IX ) IX = IX + INCX 110 CONTINUE IF( NOUNIT )& TEMP = TEMP/A( KPLUS1, J ) X( JX ) = TEMP JX = JX + INCX IF( J.GT.K )& KX = KX + INCX 120 CONTINUE END IF ELSE IF( INCX.EQ.1 )THEN DO 140, J = N, 1, -1 TEMP = X( J ) L = 1 - J DO 130, I = MIN( N, J + K ), J + 1, -1 TEMP = TEMP - A( L + I, J )*X( I ) 130 CONTINUE IF( NOUNIT )& TEMP = TEMP/A( 1, J ) X( J ) = TEMP 140 CONTINUE ELSE KX = KX + ( N - 1 )*INCX JX = KX DO 160, J = N, 1, -1 TEMP = X( JX ) IX = KX L = 1 - J DO 150, I = MIN( N, J + K ), J + 1, -1 TEMP = TEMP - A( L + I, J )*X( IX ) IX = IX - INCX 150 CONTINUE IF( NOUNIT )& TEMP = TEMP/A( 1, J ) X( JX ) = TEMP JX = JX - INCX IF( ( N - J ).GE.K )& KX = KX - INCX 160 CONTINUE END IF END IF END IF ! RETURN ! ! End of DTBSV . ! END subroutine dscal(n,da,dx,incx) ! ! scales a vector by a constant. ! uses unrolled loops for increment equal to one. ! jack dongarra, linpack, 3/11/78. ! modified to correct problem with negative increment, 8/21/90. ! double precision da,dx(1) integer i,incx,ix,m,mp1,n ! if(n.le.0)return if(incx.eq.1)go to 20 ! ! code for increment not equal to 1 ! ix = 1 if(incx.lt.0)ix = (-n+1)*incx + 1 do 10 i = 1,n dx(ix) = da*dx(ix) ix = ix + incx 10 continue return ! ! code for increment equal to 1 !c ! ! clean-up loop ! 20 m = mod(n,5) if( m .eq. 0 ) go to 40 do 30 i = 1,m dx(i) = da*dx(i) 30 continue if( n .lt. 5 ) return 40 mp1 = m + 1 do 50 i = mp1,n,5 dx(i) = da*dx(i) dx(i + 1) = da*dx(i + 1) dx(i + 2) = da*dx(i + 2) dx(i + 3) = da*dx(i + 3) dx(i + 4) = da*dx(i + 4) 50 continue return end SUBROUTINE DSYR ( UPLO, N, ALPHA, X, INCX, A, LDA ) ! .. Scalar Arguments .. DOUBLE PRECISION ALPHA INTEGER INCX, LDA, N CHARACTER*1 UPLO ! .. Array Arguments .. DOUBLE PRECISION A( LDA, * ), X( * ) ! .. ! ! Purpose ! ======= ! ! DSYR performs the symmetric rank 1 operation ! ! A := alpha*x*x' + A, ! ! where alpha is a real scalar, x is an n element vector and A is an ! n by n symmetric matrix. ! ! Parameters ! ========== ! ! UPLO - CHARACTER*1. ! On entry, UPLO specifies whether the upper or lower ! triangular part of the array A is to be referenced as ! follows: ! ! UPLO = 'U' or 'u' Only the upper triangular part of A ! is to be referenced. ! ! UPLO = 'L' or 'l' Only the lower triangular part of A ! is to be referenced. ! ! Unchanged on exit. ! ! N - INTEGER. ! On entry, N specifies the order of the matrix A. ! N must be at least zero. ! Unchanged on exit. ! ! ALPHA - DOUBLE PRECISION. ! On entry, ALPHA specifies the scalar alpha. ! Unchanged on exit. ! ! X - DOUBLE PRECISION array of dimension at least ! ( 1 + ( n - 1 )*abs( INCX ) ). ! Before entry, the incremented array X must contain the n ! element vector x. ! Unchanged on exit. ! ! INCX - INTEGER. ! On entry, INCX specifies the increment for the elements of ! X. INCX must not be zero. ! Unchanged on exit. ! ! A - DOUBLE PRECISION array of DIMENSION ( LDA, n ). ! Before entry with UPLO = 'U' or 'u', the leading n by n ! upper triangular part of the array A must contain the upper ! triangular part of the symmetric matrix and the strictly ! lower triangular part of A is not referenced. On exit, the ! upper triangular part of the array A is overwritten by the ! upper triangular part of the updated matrix. ! Before entry with UPLO = 'L' or 'l', the leading n by n ! lower triangular part of the array A must contain the lower ! triangular part of the symmetric matrix and the strictly ! upper triangular part of A is not referenced. On exit, the ! lower triangular part of the array A is overwritten by the ! lower triangular part of the updated matrix. ! ! LDA - INTEGER. ! On entry, LDA specifies the first dimension of A as declared ! in the calling (sub) program. LDA must be at least ! max( 1, n ). ! Unchanged on exit. ! ! ! Level 2 Blas routine. ! ! -- Written on 22-October-1986. ! Jack Dongarra, Argonne National Lab. ! Jeremy Du Croz, Nag Central Office. ! Sven Hammarling, Nag Central Office. ! Richard Hanson, Sandia National Labs. ! ! ! .. Parameters .. DOUBLE PRECISION ZERO PARAMETER ( ZERO = 0.0D+0 ) ! .. Local Scalars .. DOUBLE PRECISION TEMP INTEGER I, INFO, IX, J, JX, KX ! .. External Functions .. LOGICAL LSAME EXTERNAL LSAME ! .. External Subroutines .. EXTERNAL XERBLA ! .. Intrinsic Functions .. INTRINSIC MAX ! .. ! .. Executable Statements .. ! ! Test the input parameters. ! INFO = 0 IF ( .NOT.LSAME( UPLO, 'U' ).AND.& .NOT.LSAME( UPLO, 'L' ) )THEN INFO = 1 ELSE IF( N.LT.0 )THEN INFO = 2 ELSE IF( INCX.EQ.0 )THEN INFO = 5 ELSE IF( LDA.LT.MAX( 1, N ) )THEN INFO = 7 END IF IF( INFO.NE.0 )THEN CALL XERBLA( 'DSYR ', INFO ) RETURN END IF ! ! Quick return if possible. ! IF( ( N.EQ.0 ).OR.( ALPHA.EQ.ZERO ) )& RETURN ! ! Set the start point in X if the increment is not unity. ! IF( INCX.LE.0 )THEN KX = 1 - ( N - 1 )*INCX ELSE IF( INCX.NE.1 )THEN KX = 1 END IF ! ! Start the operations. In this version the elements of A are ! accessed sequentially with one pass through the triangular part ! of A. ! IF( LSAME( UPLO, 'U' ) )THEN ! ! Form A when A is stored in upper triangle. ! IF( INCX.EQ.1 )THEN DO 20, J = 1, N IF( X( J ).NE.ZERO )THEN TEMP = ALPHA*X( J ) DO 10, I = 1, J A( I, J ) = A( I, J ) + X( I )*TEMP 10 CONTINUE END IF 20 CONTINUE ELSE JX = KX DO 40, J = 1, N IF( X( JX ).NE.ZERO )THEN TEMP = ALPHA*X( JX ) IX = KX DO 30, I = 1, J A( I, J ) = A( I, J ) + X( IX )*TEMP IX = IX + INCX 30 CONTINUE END IF JX = JX + INCX 40 CONTINUE END IF ELSE ! ! Form A when A is stored in lower triangle. ! IF( INCX.EQ.1 )THEN DO 60, J = 1, N IF( X( J ).NE.ZERO )THEN TEMP = ALPHA*X( J ) DO 50, I = J, N A( I, J ) = A( I, J ) + X( I )*TEMP 50 CONTINUE END IF 60 CONTINUE ELSE JX = KX DO 80, J = 1, N IF( X( JX ).NE.ZERO )THEN TEMP = ALPHA*X( JX ) IX = JX DO 70, I = J, N A( I, J ) = A( I, J ) + X( IX )*TEMP IX = IX + INCX 70 CONTINUE END IF JX = JX + INCX 80 CONTINUE END IF END IF ! RETURN ! ! End of DSYR . ! END double precision function ddot(n,dx,incx,dy,incy) ! ! forms the dot product of two vectors. ! uses unrolled loops for increments equal to one. ! jack dongarra, linpack, 3/11/78. ! double precision dx(1),dy(1),dtemp integer i,incx,incy,ix,iy,m,mp1,n ! ddot = 0.0d0 dtemp = 0.0d0 if(n.le.0)return if(incx.eq.1.and.incy.eq.1)go to 20 ! ! code for unequal increments or equal increments ! not equal to 1 ! ix = 1 iy = 1 if(incx.lt.0)ix = (-n+1)*incx + 1 if(incy.lt.0)iy = (-n+1)*incy + 1 do 10 i = 1,n dtemp = dtemp + dx(ix)*dy(iy) ix = ix + incx iy = iy + incy 10 continue ddot = dtemp return ! ! code for both increments equal to 1 ! ! ! clean-up loop ! 20 m = mod(n,5) if( m .eq. 0 ) go to 40 do 30 i = 1,m dtemp = dtemp + dx(i)*dy(i) 30 continue if( n .lt. 5 ) go to 60 40 mp1 = m + 1 do 50 i = mp1,n,5 dtemp = dtemp + dx(i)*dy(i) + dx(i + 1)*dy(i + 1) + & dx(i + 2)*dy(i + 2) + dx(i + 3)*dy(i + 3) + dx(i + 4)*dy(i + 4) 50 continue 60 ddot = dtemp return end SUBROUTINE DGEMV ( TRANS, M, N, ALPHA, A, LDA, X, INCX,& BETA, Y, INCY ) ! .. Scalar Arguments .. DOUBLE PRECISION ALPHA, BETA INTEGER INCX, INCY, LDA, M, N CHARACTER*1 TRANS ! .. Array Arguments .. DOUBLE PRECISION A( LDA, * ), X( * ), Y( * ) ! .. ! ! Purpose ! ======= ! ! DGEMV performs one of the matrix-vector operations ! ! y := alpha*A*x + beta*y, or y := alpha*A'*x + beta*y, ! ! where alpha and beta are scalars, x and y are vectors and A is an ! m by n matrix. ! ! Parameters ! ========== ! ! TRANS - CHARACTER*1. ! On entry, TRANS specifies the operation to be performed as ! follows: ! ! TRANS = 'N' or 'n' y := alpha*A*x + beta*y. ! ! TRANS = 'T' or 't' y := alpha*A'*x + beta*y. ! ! TRANS = 'C' or 'c' y := alpha*A'*x + beta*y. ! ! Unchanged on exit. ! ! M - INTEGER. ! On entry, M specifies the number of rows of the matrix A. ! M must be at least zero. ! Unchanged on exit. ! ! N - INTEGER. ! On entry, N specifies the number of columns of the matrix A. ! N must be at least zero. ! Unchanged on exit. ! ! ALPHA - DOUBLE PRECISION. ! On entry, ALPHA specifies the scalar alpha. ! Unchanged on exit. ! ! A - DOUBLE PRECISION array of DIMENSION ( LDA, n ). ! Before entry, the leading m by n part of the array A must ! contain the matrix of coefficients. ! Unchanged on exit. ! ! LDA - INTEGER. ! On entry, LDA specifies the first dimension of A as declared ! in the calling (sub) program. LDA must be at least ! max( 1, m ). ! Unchanged on exit. ! ! X - DOUBLE PRECISION array of DIMENSION at least ! ( 1 + ( n - 1 )*abs( INCX ) ) when TRANS = 'N' or 'n' ! and at least ! ( 1 + ( m - 1 )*abs( INCX ) ) otherwise. ! Before entry, the incremented array X must contain the ! vector x. ! Unchanged on exit. ! ! INCX - INTEGER. ! On entry, INCX specifies the increment for the elements of ! X. INCX must not be zero. ! Unchanged on exit. ! ! BETA - DOUBLE PRECISION. ! On entry, BETA specifies the scalar beta. When BETA is ! supplied as zero then Y need not be set on input. ! Unchanged on exit. ! ! Y - DOUBLE PRECISION array of DIMENSION at least ! ( 1 + ( m - 1 )*abs( INCY ) ) when TRANS = 'N' or 'n' ! and at least ! ( 1 + ( n - 1 )*abs( INCY ) ) otherwise. ! Before entry with BETA non-zero, the incremented array Y ! must contain the vector y. On exit, Y is overwritten by the ! updated vector y. ! ! INCY - INTEGER. ! On entry, INCY specifies the increment for the elements of ! Y. INCY must not be zero. ! Unchanged on exit. ! ! ! Level 2 Blas routine. ! ! -- Written on 22-October-1986. ! Jack Dongarra, Argonne National Lab. ! Jeremy Du Croz, Nag Central Office. ! Sven Hammarling, Nag Central Office. ! Richard Hanson, Sandia National Labs. ! ! ! .. Parameters .. DOUBLE PRECISION ONE , ZERO PARAMETER ( ONE = 1.0D+0, ZERO = 0.0D+0 ) ! .. Local Scalars .. DOUBLE PRECISION TEMP INTEGER I, INFO, IX, IY, J, JX, JY, KX, KY, LENX, LENY ! .. External Functions .. LOGICAL LSAME EXTERNAL LSAME ! .. External Subroutines .. EXTERNAL XERBLA ! .. Intrinsic Functions .. INTRINSIC MAX ! .. ! .. Executable Statements .. ! ! Test the input parameters. ! INFO = 0 IF ( .NOT.LSAME( TRANS, 'N' ).AND.& .NOT.LSAME( TRANS, 'T' ).AND.& .NOT.LSAME( TRANS, 'C' ) )THEN INFO = 1 ELSE IF( M.LT.0 )THEN INFO = 2 ELSE IF( N.LT.0 )THEN INFO = 3 ELSE IF( LDA.LT.MAX( 1, M ) )THEN INFO = 6 ELSE IF( INCX.EQ.0 )THEN INFO = 8 ELSE IF( INCY.EQ.0 )THEN INFO = 11 END IF IF( INFO.NE.0 )THEN CALL XERBLA( 'DGEMV ', INFO ) RETURN END IF ! ! Quick return if possible. ! IF( ( M.EQ.0 ).OR.( N.EQ.0 ).OR.& ( ( ALPHA.EQ.ZERO ).AND.( BETA.EQ.ONE ) ) )& RETURN ! ! Set LENX and LENY, the lengths of the vectors x and y, and set ! up the start points in X and Y. ! IF( LSAME( TRANS, 'N' ) )THEN LENX = N LENY = M ELSE LENX = M LENY = N END IF IF( INCX.GT.0 )THEN KX = 1 ELSE KX = 1 - ( LENX - 1 )*INCX END IF IF( INCY.GT.0 )THEN KY = 1 ELSE KY = 1 - ( LENY - 1 )*INCY END IF ! ! Start the operations. In this version the elements of A are ! accessed sequentially with one pass through A. ! ! First form y := beta*y. ! IF( BETA.NE.ONE )THEN IF( INCY.EQ.1 )THEN IF( BETA.EQ.ZERO )THEN DO 10, I = 1, LENY Y( I ) = ZERO 10 CONTINUE ELSE DO 20, I = 1, LENY Y( I ) = BETA*Y( I ) 20 CONTINUE END IF ELSE IY = KY IF( BETA.EQ.ZERO )THEN DO 30, I = 1, LENY Y( IY ) = ZERO IY = IY + INCY 30 CONTINUE ELSE DO 40, I = 1, LENY Y( IY ) = BETA*Y( IY ) IY = IY + INCY 40 CONTINUE END IF END IF END IF IF( ALPHA.EQ.ZERO )& RETURN IF( LSAME( TRANS, 'N' ) )THEN ! ! Form y := alpha*A*x + y. ! JX = KX IF( INCY.EQ.1 )THEN DO 60, J = 1, N IF( X( JX ).NE.ZERO )THEN TEMP = ALPHA*X( JX ) DO 50, I = 1, M Y( I ) = Y( I ) + TEMP*A( I, J ) 50 CONTINUE END IF JX = JX + INCX 60 CONTINUE ELSE DO 80, J = 1, N IF( X( JX ).NE.ZERO )THEN TEMP = ALPHA*X( JX ) IY = KY DO 70, I = 1, M Y( IY ) = Y( IY ) + TEMP*A( I, J ) IY = IY + INCY 70 CONTINUE END IF JX = JX + INCX 80 CONTINUE END IF ELSE ! ! Form y := alpha*A'*x + y. ! JY = KY IF( INCX.EQ.1 )THEN DO 100, J = 1, N TEMP = ZERO DO 90, I = 1, M TEMP = TEMP + A( I, J )*X( I ) 90 CONTINUE Y( JY ) = Y( JY ) + ALPHA*TEMP JY = JY + INCY 100 CONTINUE ELSE DO 120, J = 1, N TEMP = ZERO IX = KX DO 110, I = 1, M TEMP = TEMP + A( I, J )*X( IX ) IX = IX + INCX 110 CONTINUE Y( JY ) = Y( JY ) + ALPHA*TEMP JY = JY + INCY 120 CONTINUE END IF END IF ! RETURN ! ! End of DGEMV . ! END INTEGER FUNCTION ILAENV( ISPEC, NAME, OPTS, N1, N2, N3,& N4 ) ! ! -- LAPACK auxiliary routine (preliminary version) -- ! Univ. of Tennessee, Univ. of California Berkeley, NAG Ltd., ! Courant Institute, Argonne National Lab, and Rice University ! February 20, 1992 ! ! .. Scalar Arguments .. CHARACTER*( * ) NAME, OPTS INTEGER ISPEC, N1, N2, N3, N4 ! .. ! ! Purpose ! ======= ! ! ILAENV is called from the LAPACK routines to choose problem-dependent ! parameters for the local environment. See ISPEC for a description of ! the parameters. ! ! This version provides a set of parameters which should give good, ! but not optimal, performance on many of the currently available ! computers. Users are encouraged to modify this subroutine to set ! the tuning parameters for their particular machine using the option ! and problem size information in the arguments. ! ! This routine will not function correctly if it is converted to all ! lower case. Converting it to all upper case is allowed. ! ! Arguments ! ========= ! ! ISPEC (input) INTEGER ! Specifies the parameter to be returned as the value of ! ILAENV. ! = 1: the optimal blocksize; if this value is 1, an unblocked ! algorithm will give the best performance. ! = 2: the minimum block size for which the block routine ! should be used; if the usable block size is less than ! this value, an unblocked routine should be used. ! = 3: the crossover point (in a block routine, for N less ! than this value, an unblocked routine should be used) ! = 4: the number of shifts, used in the nonsymmetric ! eigenvalue routines ! = 5: the minimum column dimension for blocking to be used; ! rectangular blocks must have dimension at least k by m, ! where k is given by ILAENV(2,...) and m by ILAENV(5,...) ! = 6: the crossover point for the SVD (when reducing an m by n ! matrix to bidiagonal form, if max(m,n)/min(m,n) exceeds ! this value, a QR factorization is used first to reduce ! the matrix to a triangular form.) ! = 7: the number of processors ! = 8: the crossover point for the multishift QR and QZ methods ! for nonsymmetric eigenvalue problems. ! ! NAME (input) CHARACTER*(*) ! The name of the calling subroutine, in either upper case or ! lower case. ! ! OPTS (input) CHARACTER*(*) ! The character options to the subroutine NAME, concatenated ! into a single character string. For example, UPLO = 'U', ! TRANS = 'T', and DIAG = 'N' for a triangular routine would ! be specified as OPTS = 'UTN'. ! ! N1 (input) INTEGER ! N2 (input) INTEGER ! N3 (input) INTEGER ! N4 (input) INTEGER ! Problem dimensions for the subroutine NAME; these may not all ! be required. ! ! (ILAENV) (output) INTEGER ! >= 0: the value of the parameter specified by ISPEC ! < 0: if ILAENV = -k, the k-th argument had an illegal value. ! ! Further Details ! =============== ! ! The following conventions have been used when calling ILAENV from the ! LAPACK routines: ! 1) OPTS is a concatenation of all of the character options to ! subroutine NAME, in the same order that they appear in the ! argument list for NAME, even if they are not used in determining ! the value of the parameter specified by ISPEC. ! 2) The problem dimensions N1, N2, N3, N4 are specified in the order ! that they appear in the argument list for NAME. N1 is used ! first, N2 second, and so on, and unused problem dimensions are ! passed a value of -1. ! 3) The parameter value returned by ILAENV is checked for validity in ! the calling subroutine. For example, ILAENV is used to retrieve ! the optimal blocksize for STRTRI as follows: ! ! NB = ILAENV( 1, 'STRTRI', UPLO // DIAG, N, -1, -1, -1 ) ! IF( NB.LE.1 ) NB = MAX( 1, N ) ! ! ===================================================================== ! ! .. Local Scalars .. LOGICAL CNAME, SNAME CHARACTER*1 C1 CHARACTER*2 C2, C4 CHARACTER*3 C3 CHARACTER*6 SUBNAM INTEGER I, IC, IZ, NB, NBMIN, NX ! .. ! .. Intrinsic Functions .. INTRINSIC CHAR, ICHAR, INT, MIN, REAL ! .. ! .. Executable Statements .. ! GO TO ( 100, 100, 100, 400, 500, 600, 700, 800 ) ISPEC ! ! Invalid value for ISPEC ! ILAENV = -1 RETURN ! 100 CONTINUE ! ! Convert NAME to upper case if the first character is lower case. ! ILAENV = 1 SUBNAM = NAME IC = ICHAR( SUBNAM( 1:1 ) ) IZ = ICHAR( 'Z' ) IF( IZ.EQ.90 .OR. IZ.EQ.122 ) THEN ! ! ASCII character set ! IF( IC.GE.97 .AND. IC.LE.122 ) THEN SUBNAM( 1:1 ) = CHAR( IC-32 ) DO 10 I = 2, 6 IC = ICHAR( SUBNAM( I:I ) ) IF( IC.GE.97 .AND. IC.LE.122 )& SUBNAM( I:I ) = CHAR( IC-32 ) 10 CONTINUE END IF ! ELSE IF( IZ.EQ.233 .OR. IZ.EQ.169 ) THEN ! ! EBCDIC character set ! IF( ( IC.GE.129 .AND. IC.LE.137 ) .OR.& ( IC.GE.145 .AND. IC.LE.153 ) .OR.& ( IC.GE.162 .AND. IC.LE.169 ) ) THEN SUBNAM( 1:1 ) = CHAR( IC+64 ) DO 20 I = 2, 6 IC = ICHAR( SUBNAM( I:I ) ) IF( ( IC.GE.129 .AND. IC.LE.137 ) .OR.& ( IC.GE.145 .AND. IC.LE.153 ) .OR.& ( IC.GE.162 .AND. IC.LE.169 ) )& SUBNAM( I:I ) = CHAR( IC+64 ) 20 CONTINUE END IF ! ELSE IF( IZ.EQ.218 .OR. IZ.EQ.250 ) THEN ! ! Prime machines: ASCII+128 ! IF( IC.GE.225 .AND. IC.LE.250 ) THEN SUBNAM( 1:1 ) = CHAR( IC-32 ) DO 30 I = 2, 6 IC = ICHAR( SUBNAM( I:I ) ) IF( IC.GE.225 .AND. IC.LE.250 )& SUBNAM( I:I ) = CHAR( IC-32 ) 30 CONTINUE END IF END IF ! C1 = SUBNAM( 1:1 ) SNAME = C1.EQ.'S' .OR. C1.EQ.'D' CNAME = C1.EQ.'C' .OR. C1.EQ.'Z' IF( .NOT.( CNAME .OR. SNAME ) )& RETURN C2 = SUBNAM( 2:3 ) C3 = SUBNAM( 4:6 ) C4 = C3( 2:3 ) ! GO TO ( 110, 200, 300 ) ISPEC ! 110 CONTINUE ! ! ISPEC = 1: block size ! ! In these examples, separate code is provided for setting NB for ! real and complex. We assume that NB will take the same value in ! single or double precision. ! NB = 1 ! IF( C2.EQ.'GE' ) THEN IF( C3.EQ.'TRF' ) THEN IF( SNAME ) THEN NB = 64 ELSE NB = 64 END IF ELSE IF( C3.EQ.'QRF' .OR. C3.EQ.'RQF' .OR. C3.EQ.'LQF' .OR.& C3.EQ.'QLF' ) THEN IF( SNAME ) THEN NB = 32 ELSE NB = 32 END IF ELSE IF( C3.EQ.'HRD' ) THEN IF( SNAME ) THEN NB = 32 ELSE NB = 32 END IF ELSE IF( C3.EQ.'BRD' ) THEN IF( SNAME ) THEN NB = 32 ELSE NB = 32 END IF ELSE IF( C3.EQ.'TRI' ) THEN IF( SNAME ) THEN NB = 64 ELSE NB = 64 END IF END IF ELSE IF( C2.EQ.'PO' ) THEN IF( C3.EQ.'TRF' ) THEN IF( SNAME ) THEN NB = 64 ELSE NB = 64 END IF END IF ELSE IF( C2.EQ.'SY' ) THEN IF( C3.EQ.'TRF' ) THEN IF( SNAME ) THEN NB = 64 ELSE NB = 64 END IF ELSE IF( SNAME .AND. C3.EQ.'TRD' ) THEN NB = 1 ELSE IF( SNAME .AND. C3.EQ.'GST' ) THEN NB = 64 END IF ELSE IF( CNAME .AND. C2.EQ.'HE' ) THEN IF( C3.EQ.'TRF' ) THEN NB = 64 ELSE IF( C3.EQ.'TRD' ) THEN NB = 1 ELSE IF( C3.EQ.'GST' ) THEN NB = 64 END IF ELSE IF( SNAME .AND. C2.EQ.'OR' ) THEN IF( C3( 1:1 ).EQ.'G' ) THEN IF( C4.EQ.'QR' .OR. C4.EQ.'RQ' .OR. C4.EQ.'LQ' .OR.& C4.EQ.'QL' .OR. C4.EQ.'HR' .OR. C4.EQ.'TR' .OR.& C4.EQ.'BR' ) THEN NB = 32 END IF ELSE IF( C3( 1:1 ).EQ.'M' ) THEN IF( C4.EQ.'QR' .OR. C4.EQ.'RQ' .OR. C4.EQ.'LQ' .OR.& C4.EQ.'QL' .OR. C4.EQ.'HR' .OR. C4.EQ.'TR' .OR.& C4.EQ.'BR' ) THEN NB = 32 END IF END IF ELSE IF( CNAME .AND. C2.EQ.'UN' ) THEN IF( C3( 1:1 ).EQ.'G' ) THEN IF( C4.EQ.'QR' .OR. C4.EQ.'RQ' .OR. C4.EQ.'LQ' .OR.& C4.EQ.'QL' .OR. C4.EQ.'HR' .OR. C4.EQ.'TR' .OR.& C4.EQ.'BR' ) THEN NB = 32 END IF ELSE IF( C3( 1:1 ).EQ.'M' ) THEN IF( C4.EQ.'QR' .OR. C4.EQ.'RQ' .OR. C4.EQ.'LQ' .OR.& C4.EQ.'QL' .OR. C4.EQ.'HR' .OR. C4.EQ.'TR' .OR.& C4.EQ.'BR' ) THEN NB = 32 END IF END IF ELSE IF( C2.EQ.'GB' ) THEN IF( C3.EQ.'TRF' ) THEN IF( SNAME ) THEN IF( N4.LE.64 ) THEN NB = 1 ELSE NB = 32 END IF ELSE IF( N4.LE.64 ) THEN NB = 1 ELSE NB = 32 END IF END IF END IF ELSE IF( C2.EQ.'PB' ) THEN IF( C3.EQ.'TRF' ) THEN IF( SNAME ) THEN IF( N2.LE.64 ) THEN NB = 1 ELSE NB = 32 END IF ELSE IF( N2.LE.64 ) THEN NB = 1 ELSE NB = 32 END IF END IF END IF ELSE IF( C2.EQ.'TR' ) THEN IF( C3.EQ.'TRI' ) THEN IF( SNAME ) THEN NB = 64 ELSE NB = 64 END IF END IF ELSE IF( C2.EQ.'LA' ) THEN IF( C3.EQ.'UUM' ) THEN IF( SNAME ) THEN NB = 64 ELSE NB = 64 END IF END IF ELSE IF( SNAME .AND. C2.EQ.'ST' ) THEN IF( C3.EQ.'EBZ' ) THEN NB = 1 END IF END IF ILAENV = NB RETURN ! 200 CONTINUE ! ! ISPEC = 2: minimum block size ! NBMIN = 2 IF( C2.EQ.'GE' ) THEN IF( C3.EQ.'QRF' .OR. C3.EQ.'RQF' .OR. C3.EQ.'LQF' .OR.& C3.EQ.'QLF' ) THEN IF( SNAME ) THEN NBMIN = 2 ELSE NBMIN = 2 END IF ELSE IF( C3.EQ.'HRD' ) THEN IF( SNAME ) THEN NBMIN = 2 ELSE NBMIN = 2 END IF ELSE IF( C3.EQ.'BRD' ) THEN IF( SNAME ) THEN NBMIN = 2 ELSE NBMIN = 2 END IF ELSE IF( C3.EQ.'TRI' ) THEN IF( SNAME ) THEN NBMIN = 2 ELSE NBMIN = 2 END IF END IF ELSE IF( C2.EQ.'SY' ) THEN IF( C3.EQ.'TRF' ) THEN IF( SNAME ) THEN NBMIN = 2 ELSE NBMIN = 2 END IF ELSE IF( SNAME .AND. C3.EQ.'TRD' ) THEN NBMIN = 2 END IF ELSE IF( CNAME .AND. C2.EQ.'HE' ) THEN IF( C3.EQ.'TRD' ) THEN NBMIN = 2 END IF ELSE IF( SNAME .AND. C2.EQ.'OR' ) THEN IF( C3( 1:1 ).EQ.'G' ) THEN IF( C4.EQ.'QR' .OR. C4.EQ.'RQ' .OR. C4.EQ.'LQ' .OR.& C4.EQ.'QL' .OR. C4.EQ.'HR' .OR. C4.EQ.'TR' .OR.& C4.EQ.'BR' ) THEN NBMIN = 2 END IF ELSE IF( C3( 1:1 ).EQ.'M' ) THEN IF( C4.EQ.'QR' .OR. C4.EQ.'RQ' .OR. C4.EQ.'LQ' .OR.& C4.EQ.'QL' .OR. C4.EQ.'HR' .OR. C4.EQ.'TR' .OR.& C4.EQ.'BR' ) THEN NBMIN = 2 END IF END IF ELSE IF( CNAME .AND. C2.EQ.'UN' ) THEN IF( C3( 1:1 ).EQ.'G' ) THEN IF( C4.EQ.'QR' .OR. C4.EQ.'RQ' .OR. C4.EQ.'LQ' .OR.& C4.EQ.'QL' .OR. C4.EQ.'HR' .OR. C4.EQ.'TR' .OR.& C4.EQ.'BR' ) THEN NBMIN = 2 END IF ELSE IF( C3( 1:1 ).EQ.'M' ) THEN IF( C4.EQ.'QR' .OR. C4.EQ.'RQ' .OR. C4.EQ.'LQ' .OR.& C4.EQ.'QL' .OR. C4.EQ.'HR' .OR. C4.EQ.'TR' .OR.& C4.EQ.'BR' ) THEN NBMIN = 2 END IF END IF END IF ILAENV = NBMIN RETURN ! 300 CONTINUE ! ! ISPEC = 3: crossover point ! NX = 0 IF( C2.EQ.'GE' ) THEN IF( C3.EQ.'QRF' .OR. C3.EQ.'RQF' .OR. C3.EQ.'LQF' .OR.& C3.EQ.'QLF' ) THEN IF( SNAME ) THEN NX = 128 ELSE NX = 128 END IF ELSE IF( C3.EQ.'HRD' ) THEN IF( SNAME ) THEN NX = 128 ELSE NX = 128 END IF ELSE IF( C3.EQ.'BRD' ) THEN IF( SNAME ) THEN NX = 128 ELSE NX = 128 END IF END IF ELSE IF( C2.EQ.'SY' ) THEN IF( SNAME .AND. C3.EQ.'TRD' ) THEN NX = 1 END IF ELSE IF( CNAME .AND. C2.EQ.'HE' ) THEN IF( C3.EQ.'TRD' ) THEN NX = 1 END IF ELSE IF( SNAME .AND. C2.EQ.'OR' ) THEN IF( C3( 1:1 ).EQ.'G' ) THEN IF( C4.EQ.'QR' .OR. C4.EQ.'RQ' .OR. C4.EQ.'LQ' .OR.& C4.EQ.'QL' .OR. C4.EQ.'HR' .OR. C4.EQ.'TR' .OR.& C4.EQ.'BR' ) THEN NX = 128 END IF END IF ELSE IF( CNAME .AND. C2.EQ.'UN' ) THEN IF( C3( 1:1 ).EQ.'G' ) THEN IF( C4.EQ.'QR' .OR. C4.EQ.'RQ' .OR. C4.EQ.'LQ' .OR.& C4.EQ.'QL' .OR. C4.EQ.'HR' .OR. C4.EQ.'TR' .OR.& C4.EQ.'BR' ) THEN NX = 128 END IF END IF END IF ILAENV = NX RETURN ! 400 CONTINUE ! ! ISPEC = 4: number of shifts (used by xHSEQR) ! ILAENV = 6 RETURN ! 500 CONTINUE ! ! ISPEC = 5: minimum column dimension (not used) ! ILAENV = 2 RETURN ! 600 CONTINUE ! ! ISPEC = 6: crossover point for SVD (used by xGELSS and xGESVD) ! ILAENV = INT( REAL( MIN( N1, N2 ) )*1.6E0 ) RETURN ! 700 CONTINUE ! ! ISPEC = 7: number of processors (not used) ! ILAENV = 1 RETURN ! 800 CONTINUE ! ! ISPEC = 8: crossover point for multishift (used by xHSEQR) ! ILAENV = 50 RETURN ! ! End of ILAENV ! END SUBROUTINE XERBLA( SRNAME, INFO ) ! ! -- LAPACK auxiliary routine (version 1.1) -- ! Univ. of Tennessee, Univ. of California Berkeley, NAG Ltd., ! Courant Institute, Argonne National Lab, and Rice University ! February 29, 1992 ! ! .. Scalar Arguments .. CHARACTER*6 SRNAME INTEGER INFO ! .. ! ! Purpose ! ======= ! ! XERBLA is an error handler for the LAPACK routines. ! It is called by an LAPACK routine if an input parameter has an ! invalid value. A message is printed and execution stops. ! ! Installers may consider modifying the STOP statement in order to ! call system-specific exception-handling facilities. ! ! Arguments ! ========= ! ! SRNAME (input) CHARACTER*6 ! The name of the routine which called XERBLA. ! ! INFO (input) INTEGER ! The position of the invalid parameter in the parameter list ! of the calling routine. ! ! .. Executable Statements .. ! WRITE( *, FMT = 9999 )SRNAME, INFO ! STOP ! 9999 FORMAT( ' ** On entry to ', A6, ' parameter number ', I2, ' had ',& 'an illegal value' ) ! ! End of XERBLA ! END LOGICAL FUNCTION LSAME( CA, CB ) ! ! -- LAPACK auxiliary routine (version 1.1) -- ! Univ. of Tennessee, Univ. of California Berkeley, NAG Ltd., ! Courant Institute, Argonne National Lab, and Rice University ! February 29, 1992 ! ! .. Scalar Arguments .. CHARACTER CA, CB ! .. ! ! Purpose ! ======= ! ! LSAME returns .TRUE. if CA is the same letter as CB regardless of ! case. ! ! Arguments ! ========= ! ! CA (input) CHARACTER*1 ! CB (input) CHARACTER*1 ! CA and CB specify the single characters to be compared. ! ! .. Intrinsic Functions .. INTRINSIC ICHAR ! .. ! .. Local Scalars .. INTEGER INTA, INTB, ZCODE ! .. ! .. Executable Statements .. ! ! Test if the characters are equal ! LSAME = CA.EQ.CB IF( LSAME )& RETURN ! ! Now test for equivalence if both characters are alphabetic. ! ZCODE = ICHAR( 'Z' ) ! ! Use 'Z' rather than 'A' so that ASCII can be detected on Prime ! machines, on which ICHAR returns a value with bit 8 set. ! ICHAR('A') on Prime machines returns 193 which is the same as ! ICHAR('A') on an EBCDIC machine. ! INTA = ICHAR( CA ) INTB = ICHAR( CB ) ! IF( ZCODE.EQ.90 .OR. ZCODE.EQ.122 ) THEN ! ! ASCII is assumed - ZCODE is the ASCII code of either lower or ! upper case 'Z'. ! IF( INTA.GE.97 .AND. INTA.LE.122 ) INTA = INTA - 32 IF( INTB.GE.97 .AND. INTB.LE.122 ) INTB = INTB - 32 ! ELSE IF( ZCODE.EQ.233 .OR. ZCODE.EQ.169 ) THEN ! ! EBCDIC is assumed - ZCODE is the EBCDIC code of either lower or ! upper case 'Z'. ! IF( INTA.GE.129 .AND. INTA.LE.137 .OR.& INTA.GE.145 .AND. INTA.LE.153 .OR.& INTA.GE.162 .AND. INTA.LE.169 ) INTA = INTA + 64 IF( INTB.GE.129 .AND. INTB.LE.137 .OR.& INTB.GE.145 .AND. INTB.LE.153 .OR.& INTB.GE.162 .AND. INTB.LE.169 ) INTB = INTB + 64 ! ELSE IF( ZCODE.EQ.218 .OR. ZCODE.EQ.250 ) THEN ! ! ASCII is assumed, on Prime machines - ZCODE is the ASCII code ! plus 128 of either lower or upper case 'Z'. ! IF( INTA.GE.225 .AND. INTA.LE.250 ) INTA = INTA - 32 IF( INTB.GE.225 .AND. INTB.LE.250 ) INTB = INTB - 32 END IF LSAME = INTA.EQ.INTB ! ! RETURN ! ! End of LSAME ! END