- Solve [A]{X}={B} for {X}, then multiply VCX={C}T{X},
- Solve [A]T{Y}={C} for {Y}, then multiply VYB={Y}T{B}
The matrix [A] is nonsymmetric. The data may be real or complex.
The following input may be specified:
-
FMS Module
RN = Real Nonsymmetric
CN = Complex Nonsymmetric
- Number of equations
- Number of solution vectors
- Any FMS Parameter
C E X A M P L E 18 C C Program name: CHARACTER*10 MYNAME PARAMETER (MYNAME='EXAMPLE_18') C C Number of vectors to reduce during factoring: PARAMETER (NUMRED = 0) C C Skip operations during solving (no): PARAMETER (ISKIP = 0) C C FMS matrix and vector file attributes: INTEGER LUF(25) INTEGER LUB(25), LUC(25), LUX(25), LUY(25), LUE(25) INTEGER LUS0(25) C C FMS memory management requires the following arrays: POINTER (CMD_PTR, CMD) POINTER (RMD_PTR, RMD) POINTER (IMD_PTR, IMD) COMPLEX*16 CMD(0:1) REAL*8 RMD(0:1) INTEGER IMD(0:1) C C Input parameters: LOGICAL ASK INTEGER ASK_I C C Local variables: INTEGER IDTYPE, I, LB, LC, LX, LY, LE, LDISK, NV INTEGER LENB, LENC, LENX, LENY, LENE REAL*8 EI, ERROR, ERROR2 REAL*8 R_VCX, R_VYB, RZERO, RONE COMPLEX*16 C_VCX, C_VYB, CZERO, CONE C C Profile vector for a full matrix: INTEGER LOWEQ(1) DATA LOWEQ/-1/ C C Common block to communicate with fill routines: COMMON /MYDATA/MOD, N, NRHS C C (1) Initialize FMS: CALL FMSINI CALL FMSPSH(MYNAME) CALL FMSIST ('MDATAU', 1) RZERO = 0.0D0 RONE = 1.0D0 CZERO = DCMPLX(0.0D0, 0.0D0) CONE = DCMPLX(1.0D0, 1.0D0) 1 CONTINUE WRITE (6,*) 'Enter the FMS module number (2 or 5)' WRITE (6,*) ' 2 = Real Nonsymmetric' WRITE (6,*) ' 5 = Complex Nonsymmetric' READ (5,*) MOD IF(MOD .EQ. 2) THEN IDTYPE = 1 ELSE IF (MOD .EQ. 5) THEN IDTYPE = 2 ELSE GO TO 1 END IF N = ASK_I('Enter the number of equations') NRHS = ASK_I('Enter the number of vectors') WRITE (6,*) 'You may now alter any FMS parameter.' WRITE (6,*) 'When you are finished, type the letters RETURN' CALL FMSSET CALL FMSIGT ('MEMPTR', IMD_PTR) CALL FMSIGT ('MEMPTR', RMD_PTR) CALL FMSIGT ('MEMPTR', CMD_PTR) C C (2) Open FMS files: IF(MOD.EQ.2) CALL RNDI (LOWEQ, N, 'LUF', LUF) IF(MOD.EQ.5) CALL CNDI (LOWEQ, N, 'LUF', LUF) CALL FMSOV (N, IDTYPE, NRHS, 'LUB', LUB) CALL FMSOV (N, IDTYPE, NRHS, 'LUC', LUC) CALL FMSOV (N, IDTYPE, NRHS, 'LUX', LUX) CALL FMSOV (N, IDTYPE, NRHS, 'LUY', LUY) CALL FMSOV (N, IDTYPE, NRHS, 'LUE', LUE) C C Allocate storage for one {B}, {C}, {X}, {Y} and {E}: C Initialize {B} and {C}: IF(IDTYPE .EQ. 1) THEN LENB = LUB(4) CALL FMSRMG (RMD, LB, LENB) LENC = LUC(4) CALL FMSRMG (RMD, LC, LENC) LENX = LUX(4) CALL FMSRMG (RMD, LX, LENX) LENY = LUY(4) CALL FMSRMG (RMD, LY, LENY) LENE = LUE(4) CALL FMSRMG (RMD, LE, LENE) DO I = 0,(N-1) RMD(LB+I) = RONE RMD(LC+I) = RONE END DO ELSE LENB = LUB(4)/2 CALL FMSCMG (CMD, LB, LENB) LENC = LUC(4)/2 CALL FMSCMG (CMD, LC, LENC) LENX = LUX(4)/2 CALL FMSCMG (CMD, LX, LENX) LENY = LUY(4)/2 CALL FMSCMG (CMD, LY, LENY) LENE = LUE(4)/2 CALL FMSCMG (CMD, LE, LENE) DO I = 0,(N-1) CMD(LB+I) = CONE CMD(LC+I) = CONE END DO END IF C C (3) Write data to FMS files: C C Matrix file: IF(MOD.EQ.2) CALL RNDA (LUS0, 0, LUF) IF(MOD.EQ.5) CALL CNDA (LUS0, 0, LUF) C C Vector files: LDISK = 1 DO NV = 1,NRHS IF(IDTYPE .EQ. 1) THEN CALL FMSWRT (LUB(1), LDISK, RMD(LB), LUB(4)) CALL FMSWRT (LUC(1), LDISK, RMD(LC), LUC(4)) ELSE CALL FMSWRT (LUB(1), LDISK, CMD(LB), LUB(4)) CALL FMSWRT (LUC(1), LDISK, CMD(LC), LUC(4)) END IF LDISK = LDISK + LUB(4) END DO C C (4) Perform matrix algebra: C C Populate and factor matrix [A]: IF(MOD.EQ.2) CALL RNDF (LUF, LUF, LUB, LUX, NUMRED) IF(MOD.EQ.5) CALL CNDF (LUF, LUF, LUB, LUX, NUMRED) C C Solve [A]{X} = {B}: IF(MOD.EQ.2) CALL RNDS (LUF, LUB, LUX, NRHS, ISKIP) IF(MOD.EQ.5) CALL CNDS (LUF, LUB, LUX, NRHS, ISKIP) C C Solve [A]T{Y} = {C}: C Set the FMS Parameter LUTRAN to solve with [A]T: CALL FMSIST ('LUTRAN', 1) IF(MOD.EQ.2) CALL RNDS (LUF, LUC, LUY, NRHS, ISKIP) IF(MOD.EQ.5) CALL CNDS (LUF, LUC, LUY, NRHS, ISKIP) C Reset the FMS Parameter LUTRAN for a normal solve: CALL FMSIST ('LUTRAN', 0) C C Repopulate matrix [A]: CALL FMSIST ('MDATAU', 1) IF(MOD.EQ.2) CALL RNDA (LUS0, 0, LUF) IF(MOD.EQ.5) CALL CNDA (LUS0, 0, LUF) C C Multiply {E} = {B} - [A]{X}: CALL FMSMM (LUE, LUB, -1, 'N', LUF, LUX) C C (5) Read data from FMS files: C Compute VCX and VYB: C Check the answer: ERROR = RZERO ERROR2 = RZERO LDISK = 1 NWREAD = IDTYPE*N DO NV = 1,NRHS IF(IDTYPE .EQ. 1) THEN C Read in {B}nv, {C}nv, {X}nv, {Y}nv and {E}nv CALL FMSRED (LUB(1), LDISK, RMD(LB), NWREAD) CALL FMSRED (LUC(1), LDISK, RMD(LC), NWREAD) CALL FMSRED (LUX(1), LDISK, RMD(LX), NWREAD) CALL FMSRED (LUY(1), LDISK, RMD(LY), NWREAD) CALL FMSRED (LUE(1), LDISK, RMD(LE), NWREAD) LDISK = LDISK + LUB(4) R_VCX = RZERO R_VYB = RZERO DO I = 0,(N-1) R_VCX = R_VCX + RMD(LC+I) * RMD(LX+I) R_VYB = R_VYB + RMD(LY+I) * RMD(LB+I) EI = ABS(RMD(LE+I)) IF(EI .GT. ERROR) ERROR = EI END DO EI = ABS(R_VCX - R_VYB) IF(EI .GT. ERROR2) ERROR2 = EI ELSE C Read in {B}nv, {C}nv, {X}nv, {Y}nv and {E}nv CALL FMSRED (LUB(1), LDISK, CMD(LB), NWREAD) CALL FMSRED (LUC(1), LDISK, CMD(LC), NWREAD) CALL FMSRED (LUX(1), LDISK, CMD(LX), NWREAD) CALL FMSRED (LUY(1), LDISK, CMD(LY), NWREAD) CALL FMSRED (LUE(1), LDISK, CMD(LE), NWREAD) LDISK = LDISK + LUB(4) C_VCX = CZERO C_VYB = CZERO DO I = 0,(N-1) C_VCX = C_VCX + CMD(LC+I) * CMD(LX+I) C_VYB = C_VYB + CMD(LY+I) * CMD(LB+I) EI = ABS(CMD(LE+I)) IF(EI .GT. ERROR) ERROR = EI END DO EI = ABS(C_VCX - C_VYB) IF(EI .GT. ERROR2) ERROR2 = EI END IF END DO WRITE(6,*) 'MAXIMUM ERROR IN {B}-[A]{X} =', ERROR WRITE(6,*) 'MAXIMUM ERROR IN {C}T{X}-{Y}T{B} =', ERROR2 C C (6) Close FMS files: CALL FMSCV (LUE) CALL FMSCV (LUY) CALL FMSCV (LUX) CALL FMSCV (LUC) CALL FMSCV (LUB) CALL FMSCM (LUF) C C Return the memory: IF(IDTYPE .EQ. 1) THEN CALL FMSRMR (RMD, LE, LENE) CALL FMSRMR (RMD, LY, LENY) CALL FMSRMR (RMD, LX, LENX) CALL FMSRMR (RMD, LC, LENC) CALL FMSRMR (RMD, LB, LENB) ELSE CALL FMSCMR (CMD, LE, LENE) CALL FMSCMR (CMD, LY, LENY) CALL FMSCMR (CMD, LX, LENX) CALL FMSCMR (CMD, LC, LENC) CALL FMSCMR (CMD, LB, LENB) END IF IF(ASK('Do you want another solution?')) GO TO 1 CALL FMSPOP(MYNAME) CALL FMSEND END C======================================================================= LOGICAL FUNCTION ASK(QUESTION) C======================================================================= CHARACTER* (*) QUESTION CHARACTER*1 IYN WRITE(6,2000) QUESTION READ (5,1000) IYN IF( (IYN .EQ. 'Y') .OR. (IYN .EQ. 'y') ) THEN ASK = .TRUE. ELSE ASK = .FALSE. END IF RETURN 1000 FORMAT (A) 2000 FORMAT (1X,A,' (y,n)>') END C======================================================================= INTEGER FUNCTION ASK_I(STRING) C======================================================================= CHARACTER* (*) STRING WRITE(6,2000) STRING READ (5,*) ASK_I RETURN 2000 FORMAT (1X,A,'>') END C======================================================================= SUBROUTINE RSUBLK (A, D, LOWEQ, LOCEQ, IROW1, IROW2, 1 JCOL1, JCOL2, IJSTEP) C======================================================================= INTEGER IROW1, IROW2, JCOL1, JCOL2, IJSTEP INTEGER LOWEQ(*), LOCEQ(*) REAL*8 A(0:*), D(*), AIJ LOGICAL IDIAG C IF(IROW1 .EQ. JCOL1) THEN IDIAG = .TRUE. ELSE IDIAG = .FALSE. END IF AIJ = SQRT(2.0) DO I = IROW1,IROW2 J1 = MAX0(LOWEQ(I),JCOL1) J2 = MIN0(I-1,JCOL2) DO J = J1,J2 A(LOCEQ(I)+IJSTEP*J) = AIJ AIJ = SQRT(10.0*(AIJ - INT(AIJ))) END DO IF(IDIAG) THEN D(I) = AIJ AIJ = SQRT(10.0*(AIJ - INT(AIJ))) END IF END DO RETURN END C======================================================================= SUBROUTINE RNUBLK (A, D, LOWEQ, LOCEQ, IROW1, IROW2, 1 JCOL1, JCOL2, IJSTEP) C======================================================================= INTEGER IROW1, IROW2, JCOL1, JCOL2, IJSTEP INTEGER LOWEQ(*), LOCEQ(*) REAL*8 A(0:*), D(*), AIJ C AIJ = SQRT(2.0) DO J = JCOL1,JCOL2 I1 = MAX0(LOWEQ(J),IROW1) I2 = MIN0(J-1,IROW2) DO I = I1,I2 A(LOCEQ(J)+IJSTEP*I) = AIJ AIJ = SQRT(10.0*(AIJ - INT(AIJ))) END DO END DO RETURN END C======================================================================= SUBROUTINE CSUBLK (A, D, LOWEQ, LOCEQ, IROW1, IROW2, 1 JCOL1, JCOL2, IJSTEP) C======================================================================= INTEGER IROW1, IROW2, JCOL1, JCOL2, IJSTEP INTEGER LOWEQ(*), LOCEQ(*) REAL*8 AIJ, AIJ_RE, AIJ_IM COMPLEX*16 A(0:*), D(*) LOGICAL IDIAG C IF(IROW1 .EQ. JCOL1) THEN IDIAG = .TRUE. ELSE IDIAG = .FALSE. END IF AIJ = SQRT(2.0) DO I = IROW1,IROW2 J1 = MAX0(LOWEQ(I),JCOL1) J2 = MIN0(I-1,JCOL2) DO J = J1,J2 AIJ_RE = AIJ AIJ = SQRT(10.0*(AIJ - INT(AIJ))) AIJ_IM = AIJ A(LOCEQ(I)+IJSTEP*J) = DCMPLX(AIJ_RE,AIJ_IM) AIJ = SQRT(10.0*(AIJ - INT(AIJ))) END DO IF(IDIAG) THEN AIJ_RE = AIJ AIJ = SQRT(10.0*(AIJ - INT(AIJ))) AIJ_IM = AIJ D(I) = DCMPLX(AIJ_RE,AIJ_IM) AIJ = SQRT(10.0*(AIJ - INT(AIJ))) END IF END DO RETURN END C======================================================================= SUBROUTINE CNUBLK (A, D, LOWEQ, LOCEQ, IROW1, IROW2, 1 JCOL1, JCOL2, IJSTEP) C======================================================================= INTEGER IROW1, IROW2, JCOL1, JCOL2, IJSTEP INTEGER LOWEQ(*), LOCEQ(*) REAL*8 AIJ, AIJ_RE, AIJ_IM COMPLEX*16 A(0:*), D(*) C AIJ = SQRT(2.0) DO J = JCOL1,JCOL2 I1 = MAX0(LOWEQ(J),IROW1) I2 = MIN0(J-1,IROW2) DO I = I1,I2 AIJ_RE = AIJ AIJ = SQRT(10.0*(AIJ - INT(AIJ))) AIJ_IM = AIJ A(LOCEQ(J)+IJSTEP*I) = DCMPLX(AIJ_RE,AIJ_IM) AIJ = SQRT(10.0*(AIJ - INT(AIJ))) END DO END DO RETURN END C======================================================================= SUBROUTINE RNUSLB (A, LOCI, LOCJ, LUFLAG, JEQN1, JEQN2, NUMEQ) C======================================================================= INTEGER JEQN1, JEQN2, NUMEQ INTEGER LOCI(NUMEQ), LOCJ(JEQN1:JEQN2,2), LUFLAG(NUMEQ) REAL*8 A(0:*), AIJ C AIJ = SQRT(2.0) DO I=1,NUMEQ DO J=JEQN1,JEQN2 LIJ = LOCI(I) + LOCJ(J,LUFLAG(I)) A(LIJ) = AIJ AIJ = SQRT(10.0*(AIJ - INT(AIJ))) END DO END DO RETURN END C======================================================================= SUBROUTINE CNUSLB (A, LOCI, LOCJ, LUFLAG, JEQN1, JEQN2, NUMEQ) C======================================================================= INTEGER JEQN1, JEQN2, NUMEQ INTEGER LOCI(NUMEQ), LOCJ(JEQN1:JEQN2,2), LUFLAG(NUMEQ) COMPLEX*16 A(0:*) REAL*8 AIJ, AIJ_RE, AIJ_IM C AIJ = SQRT(2.0) DO I=1,NUMEQ DO J=JEQN1,JEQN2 LIJ = LOCI(I) + LOCJ(J,LUFLAG(I)) AIJ_RE = AIJ AIJ = SQRT(10.0*(AIJ - INT(AIJ))) AIJ_IM = AIJ A(LIJ) = DCMPLX(AIJ_RE,AIJ_IM) AIJ = SQRT(10.0*(AIJ - INT(AIJ))) END DO END DO RETURN END