*DECK SPELIP SUBROUTINE SPELIP (INTL, IORDER, A, B, M, MBDCND, BDA, ALPHA, BDB, + BETA, C, D, N, NBDCND, BDC, GAMA, BDD, XNU, COFX, COFY, AN, BN, + CN, DN, UN, ZN, AM, BM, CM, DM, UM, ZM, GRHS, USOL, IDMN, W, + PERTRB, IERROR) C***BEGIN PROLOGUE SPELIP C***SUBSIDIARY C***PURPOSE Subsidiary to SEPELI C***LIBRARY SLATEC C***TYPE SINGLE PRECISION (SPELIP-S) C***AUTHOR (UNKNOWN) C***DESCRIPTION C C SPELIP sets up vectors and arrays for input to BLKTRI C and computes a second order solution in USOL. A return jump to C SEPELI occurs if IORDER=2. If IORDER=4 a fourth order C solution is generated in USOL. C C***SEE ALSO SEPELI C***ROUTINES CALLED BLKTRI, CHKSNG, DEFER, MINSOL, ORTHOG, TRISP C***COMMON BLOCKS SPLPCM C***REVISION HISTORY (YYMMDD) C 801001 DATE WRITTEN C 890531 Changed all specific intrinsics to generic. (WRB) C 891214 Prologue converted to Version 4.0 format. (BAB) C 900402 Added TYPE section. (WRB) C***END PROLOGUE SPELIP C DIMENSION BDA(*) ,BDB(*) ,BDC(*) ,BDD(*) , 1 W(*) DIMENSION GRHS(IDMN,*) ,USOL(IDMN,*) DIMENSION AN(*) ,BN(*) ,CN(*) ,DN(*) , 1 UN(*) ,ZN(*) DIMENSION AM(*) ,BM(*) ,CM(*) ,DM(*) , 1 UM(*) ,ZM(*) COMMON /SPLPCM/ KSWX ,KSWY ,K ,L , 1 AIT ,BIT ,CIT ,DIT , 2 MIT ,NIT ,IS ,MS , 3 JS ,NS ,DLX ,DLY , 4 TDLX3 ,TDLY3 ,DLX4 ,DLY4 LOGICAL SINGLR EXTERNAL COFX ,COFY C***FIRST EXECUTABLE STATEMENT SPELIP KSWX = MBDCND+1 KSWY = NBDCND+1 K = M+1 L = N+1 AIT = A BIT = B CIT = C DIT = D C C SET RIGHT HAND SIDE VALUES FROM GRHS IN USOL ON THE INTERIOR C AND NON-SPECIFIED BOUNDARIES. C DO 20 I=2,M DO 10 J=2,N USOL(I,J) = GRHS(I,J) 10 CONTINUE 20 CONTINUE IF (KSWX.EQ.2 .OR. KSWX.EQ.3) GO TO 40 DO 30 J=2,N USOL(1,J) = GRHS(1,J) 30 CONTINUE 40 CONTINUE IF (KSWX.EQ.2 .OR. KSWX.EQ.5) GO TO 60 DO 50 J=2,N USOL(K,J) = GRHS(K,J) 50 CONTINUE 60 CONTINUE IF (KSWY.EQ.2 .OR. KSWY.EQ.3) GO TO 80 DO 70 I=2,M USOL(I,1) = GRHS(I,1) 70 CONTINUE 80 CONTINUE IF (KSWY.EQ.2 .OR. KSWY.EQ.5) GO TO 100 DO 90 I=2,M USOL(I,L) = GRHS(I,L) 90 CONTINUE 100 CONTINUE IF (KSWX.NE.2 .AND. KSWX.NE.3 .AND. KSWY.NE.2 .AND. KSWY.NE.3) 1 USOL(1,1) = GRHS(1,1) IF (KSWX.NE.2 .AND. KSWX.NE.5 .AND. KSWY.NE.2 .AND. KSWY.NE.3) 1 USOL(K,1) = GRHS(K,1) IF (KSWX.NE.2 .AND. KSWX.NE.3 .AND. KSWY.NE.2 .AND. KSWY.NE.5) 1 USOL(1,L) = GRHS(1,L) IF (KSWX.NE.2 .AND. KSWX.NE.5 .AND. KSWY.NE.2 .AND. KSWY.NE.5) 1 USOL(K,L) = GRHS(K,L) I1 = 1 C C SET SWITCHES FOR PERIODIC OR NON-PERIODIC BOUNDARIES C MP = 1 NP = 1 IF (KSWX .EQ. 1) MP = 0 IF (KSWY .EQ. 1) NP = 0 C C SET DLX,DLY AND SIZE OF BLOCK TRI-DIAGONAL SYSTEM GENERATED C IN NINT,MINT C DLX = (BIT-AIT)/M MIT = K-1 IF (KSWX .EQ. 2) MIT = K-2 IF (KSWX .EQ. 4) MIT = K DLY = (DIT-CIT)/N NIT = L-1 IF (KSWY .EQ. 2) NIT = L-2 IF (KSWY .EQ. 4) NIT = L TDLX3 = 2.0*DLX**3 DLX4 = DLX**4 TDLY3 = 2.0*DLY**3 DLY4 = DLY**4 C C SET SUBSCRIPT LIMITS FOR PORTION OF ARRAY TO INPUT TO BLKTRI C IS = 1 JS = 1 IF (KSWX.EQ.2 .OR. KSWX.EQ.3) IS = 2 IF (KSWY.EQ.2 .OR. KSWY.EQ.3) JS = 2 NS = NIT+JS-1 MS = MIT+IS-1 C C SET X - DIRECTION C DO 110 I=1,MIT XI = AIT+(IS+I-2)*DLX CALL COFX (XI,AI,BI,CI) AXI = (AI/DLX-0.5*BI)/DLX BXI = -2.*AI/DLX**2+CI CXI = (AI/DLX+0.5*BI)/DLX AM(I) = AXI BM(I) = BXI CM(I) = CXI 110 CONTINUE C C SET Y DIRECTION C DO 120 J=1,NIT YJ = CIT+(JS+J-2)*DLY CALL COFY (YJ,DJ,EJ,FJ) DYJ = (DJ/DLY-0.5*EJ)/DLY EYJ = (-2.*DJ/DLY**2+FJ) FYJ = (DJ/DLY+0.5*EJ)/DLY AN(J) = DYJ BN(J) = EYJ CN(J) = FYJ 120 CONTINUE C C ADJUST EDGES IN X DIRECTION UNLESS PERIODIC C AX1 = AM(1) CXM = CM(MIT) GO TO (170,130,150,160,140),KSWX C C DIRICHLET-DIRICHLET IN X DIRECTION C 130 AM(1) = 0.0 CM(MIT) = 0.0 GO TO 170 C C MIXED-DIRICHLET IN X DIRECTION C 140 AM(1) = 0.0 BM(1) = BM(1)+2.*ALPHA*DLX*AX1 CM(1) = CM(1)+AX1 CM(MIT) = 0.0 GO TO 170 C C DIRICHLET-MIXED IN X DIRECTION C 150 AM(1) = 0.0 AM(MIT) = AM(MIT)+CXM BM(MIT) = BM(MIT)-2.*BETA*DLX*CXM CM(MIT) = 0.0 GO TO 170 C C MIXED - MIXED IN X DIRECTION C 160 CONTINUE AM(1) = 0.0 BM(1) = BM(1)+2.*DLX*ALPHA*AX1 CM(1) = CM(1)+AX1 AM(MIT) = AM(MIT)+CXM BM(MIT) = BM(MIT)-2.*DLX*BETA*CXM CM(MIT) = 0.0 170 CONTINUE C C ADJUST IN Y DIRECTION UNLESS PERIODIC C DY1 = AN(1) FYN = CN(NIT) GO TO (220,180,200,210,190),KSWY C C DIRICHLET-DIRICHLET IN Y DIRECTION C 180 CONTINUE AN(1) = 0.0 CN(NIT) = 0.0 GO TO 220 C C MIXED-DIRICHLET IN Y DIRECTION C 190 CONTINUE AN(1) = 0.0 BN(1) = BN(1)+2.*DLY*GAMA*DY1 CN(1) = CN(1)+DY1 CN(NIT) = 0.0 GO TO 220 C C DIRICHLET-MIXED IN Y DIRECTION C 200 AN(1) = 0.0 AN(NIT) = AN(NIT)+FYN BN(NIT) = BN(NIT)-2.*DLY*XNU*FYN CN(NIT) = 0.0 GO TO 220 C C MIXED - MIXED DIRECTION IN Y DIRECTION C 210 CONTINUE AN(1) = 0.0 BN(1) = BN(1)+2.*DLY*GAMA*DY1 CN(1) = CN(1)+DY1 AN(NIT) = AN(NIT)+FYN BN(NIT) = BN(NIT)-2.0*DLY*XNU*FYN CN(NIT) = 0.0 220 IF (KSWX .EQ. 1) GO TO 270 C C ADJUST USOL ALONG X EDGE C DO 260 J=JS,NS IF (KSWX.NE.2 .AND. KSWX.NE.3) GO TO 230 USOL(IS,J) = USOL(IS,J)-AX1*USOL(1,J) GO TO 240 230 USOL(IS,J) = USOL(IS,J)+2.0*DLX*AX1*BDA(J) 240 IF (KSWX.NE.2 .AND. KSWX.NE.5) GO TO 250 USOL(MS,J) = USOL(MS,J)-CXM*USOL(K,J) GO TO 260 250 USOL(MS,J) = USOL(MS,J)-2.0*DLX*CXM*BDB(J) 260 CONTINUE 270 IF (KSWY .EQ. 1) GO TO 320 C C ADJUST USOL ALONG Y EDGE C DO 310 I=IS,MS IF (KSWY.NE.2 .AND. KSWY.NE.3) GO TO 280 USOL(I,JS) = USOL(I,JS)-DY1*USOL(I,1) GO TO 290 280 USOL(I,JS) = USOL(I,JS)+2.0*DLY*DY1*BDC(I) 290 IF (KSWY.NE.2 .AND. KSWY.NE.5) GO TO 300 USOL(I,NS) = USOL(I,NS)-FYN*USOL(I,L) GO TO 310 300 USOL(I,NS) = USOL(I,NS)-2.0*DLY*FYN*BDD(I) 310 CONTINUE 320 CONTINUE C C SAVE ADJUSTED EDGES IN GRHS IF IORDER=4 C IF (IORDER .NE. 4) GO TO 350 DO 330 J=JS,NS GRHS(IS,J) = USOL(IS,J) GRHS(MS,J) = USOL(MS,J) 330 CONTINUE DO 340 I=IS,MS GRHS(I,JS) = USOL(I,JS) GRHS(I,NS) = USOL(I,NS) 340 CONTINUE 350 CONTINUE IORD = IORDER PERTRB = 0.0 C C CHECK IF OPERATOR IS SINGULAR C CALL CHKSNG (MBDCND,NBDCND,ALPHA,BETA,GAMA,XNU,COFX,COFY,SINGLR) C C COMPUTE NON-ZERO EIGENVECTOR IN NULL SPACE OF TRANSPOSE C IF SINGULAR C IF (SINGLR) CALL TRISP (MIT,AM,BM,CM,DM,UM,ZM) IF (SINGLR) CALL TRISP (NIT,AN,BN,CN,DN,UN,ZN) C C MAKE INITIALIZATION CALL TO BLKTRI C IF (INTL .EQ. 0) 1 CALL BLKTRI (INTL,NP,NIT,AN,BN,CN,MP,MIT,AM,BM,CM,IDMN, 2 USOL(IS,JS),IERROR,W) IF (IERROR .NE. 0) RETURN C C ADJUST RIGHT HAND SIDE IF NECESSARY C 360 CONTINUE IF (SINGLR) CALL ORTHOG (USOL,IDMN,ZN,ZM,PERTRB) C C COMPUTE SOLUTION C CALL BLKTRI (I1,NP,NIT,AN,BN,CN,MP,MIT,AM,BM,CM,IDMN,USOL(IS,JS), 1 IERROR,W) IF (IERROR .NE. 0) RETURN C C SET PERIODIC BOUNDARIES IF NECESSARY C IF (KSWX .NE. 1) GO TO 380 DO 370 J=1,L USOL(K,J) = USOL(1,J) 370 CONTINUE 380 IF (KSWY .NE. 1) GO TO 400 DO 390 I=1,K USOL(I,L) = USOL(I,1) 390 CONTINUE 400 CONTINUE C C MINIMIZE SOLUTION WITH RESPECT TO WEIGHTED LEAST SQUARES C NORM IF OPERATOR IS SINGULAR C IF (SINGLR) CALL MINSOL (USOL,IDMN,ZN,ZM,PRTRB) C C RETURN IF DEFERRED CORRECTIONS AND A FOURTH ORDER SOLUTION ARE C NOT FLAGGED C IF (IORD .EQ. 2) RETURN IORD = 2 C C COMPUTE NEW RIGHT HAND SIDE FOR FOURTH ORDER SOLUTION C CALL DEFER (COFX,COFY,IDMN,USOL,GRHS) GO TO 360 END