C COPYRIGHT (C) 1983 GLENN EVERHART C PERMISSION IS GIVEN TO ANYONE TO USE, DISTRIBUTE, OR COPY THIS C PROGRAM FREELY BUT NOT TO SELL IT COMMERICALLY. C SOLVE THE MATRIX EQUATION AX=B, WHERE A AND B ARE MATRICES, C N BY N FOR A, N BY M FOR B C RETURNS VALUE OF X IN LOCATION USED FOR B ON INPUT. C THE STORAGE CONVENTION FOR FORTRAN IS C A(COL, ROW), B(COL,ROW). SUBROUTINE MTXEQU(A1,A2,B1,B2,N,M,D) C A1,A2 ARE DIMENSIONS OF A SUBMATRIX ORIGIN IN XVBLS C B1,B2 ARE DIMS OF B SUBMATRIX INCLUDE 'VKLUGPRM.FTN' C C NOTE THIS PROGRAM MUST BE MODIFIED TO WORK WITHIN THE SPREAD C SHEET MATRIX RATHER THAN JUST ASSUMING THAT THE N DIMENSION C AND M DIMENSION GIVE THE STORAGE ADDRESSES... ALTERNATIVELY, C THE PROGRAM MUST OPERATE ONLY ON COPIED, DENSELY STORED C MATRICES. C C C ORIGINAL PROGRAM TEXT FOLLOWS: C DIMENSION A(1),B(1) CC ALTER DECLARATIONS FOR USE WITH SPREAD SHEET C REAL*8 A,B C KMAX=N-1 C DO 90 K=1,KMAX C AMAX=0. C J2=K C DO 20 J1=K,N C IK=(J1-1)*N+K C IF(ABS(AMAX)-ABS(A(IK)))10,20,20 C10 AMAX=A(IK) C J2=J1 C20 CONTINUE CC EXCHANGE ROW K,J2 IF NECESSARY C IF(J2-K)30,60,30 C30 DO 40 J=K,N C J3=(K-1)*N+J C J4=(J2-1)*N+J C SAVE=A(J3) C A(J3)=A(J4) C A(J4)=SAVE C40 CONTINUE C DO 50 J=1,M C J3=(K-1)*M+J C J4=(J2-1)*M+J C SAVE=B(J3) C B(J3)=B(J4) C50 B(J4)=SAVE CC REDUCTION C60 K1=K+1 C KK=(K-1)*N+K C DO 80 I=K1,N C IK=(I-1)*N+K C DO 70 J=K1,N C IJ=(I-1)*M+J C KJ=(K-1)*M+J C70 A(IJ)=A(IJ)-A(KJ)*A(IK)/A(KK) C DO 80 J=1,M C IJ=(I-1)*M+J C KJ=(K-1)*N+J C80 B(IJ)=B(IJ)-B(KJ)*A(IK)/A(KK) C90 CONTINUE CC SUBSTITUTE BACK CC NN=(N-1)*N+N C NN=N*N C DO 110 J=1,M C NJ=(N-1)*M+J C B(NJ)=B(NJ)/A(NN) C I1MAX=N-1 C IF(I1MAX)110,110,95 C95 DO 111 I1=1,I1MAX C I=N-I1 C IJ=(I-1)*M+J C II=(I-1)*N+I C I2=I+1 C DO 100 L=I2,N C IL=(I-1)*N+L C LJ=(L-1)*M+J C100 B(IJ)=B(IJ)-A(IL)*B(LJ) C B(IJ)=B(IJ)/A(II) C111 CONTINUE C110 CONTINUE C RETURN C END INTEGER A1,A2,B1,B2 C DIMENSION A(1),B(1) C ALTER DECLARATIONS FOR USE WITH SPREAD SHEET C NOTE THAT OUR COLUMN DIMENSION IS RRW, NOT N OR M HERE C SUBSCRIPTS ARE (ROW-1)*COL-DIMENSION + COL C THEREFORE, CHANGE *N OR *M IN SUBSCRIPT COMPUTATIONS TO C *RRW REAL*8 A,B,AW1,AW2,BW1,BW2,AW3,AW4 INTEGER ABASE,BBASE ABASE=(A2-1)*RRW+A1-1 BBASE=(B2-1)*RRW+B1-1 D=1. KMAX=N-1 DO 90 K=1,KMAX AMAX=0. J2=K DO 20 J1=K,N IK=(J1-1)*RRW+K CALL XVBLGT(IK+ABASE,0,A) IF(ABS(AMAX)-ABS(A))10,20,20 10 AMAX=A J2=J1 20 CONTINUE C EXCHANGE ROW K,J2 IF NECESSARY IF(J2-K)30,60,30 30 DO 40 J=K,N J3=(K-1)*RRW+J J4=(J2-1)*RRW+J CALL XVBLGT(J3+ABASE,0,SAVE) C SAVE=A(J3) CALL XVBLGT(J4+ABASE,0,AW1) CALL XVBLST(J3+ABASE,0,AW1) CALL XVBLST(J4+ABASE,0,SAVE) C A(J3)=A(J4) C A(J4)=SAVE 40 CONTINUE DO 50 J=1,M J3=(K-1)*RRW+J J4=(J2-1)*RRW+J C SAVE=B(J3) C B(J3)=B(J4) C50 B(J4)=SAVE CALL XVBLGT(J3+BBASE,0,SAVE) CALL XVBLGT(J4+BBASE,0,BW1) CALL XVBLST(J3+BBASE,0,BW1) CALL XVBLST(J4+BBASE,0,SAVE) 50 CONTINUE C REDUCTION 60 K1=K+1 KK=(K-1)*RRW+K CALL XVBLGT(KK+ABASE,0,A) IF(A.EQ.0)GOTO 999 C IF(A(KK).EQ.0.)GOTO 999 DO 80 I=K1,N IK=(I-1)*RRW+K DO 70 J=K1,N IJ=(I-1)*RRW+J KJ=(K-1)*RRW+J C70 A(IJ)=A(IJ)-A(KJ)*A(IK)/A(KK) CALL XVBLGT(IJ+ABASE,0,AW1) CALL XVBLGT(KJ+ABASE,0,AW2) CALL XVBLGT(IK+ABASE,0,AW3) CALL XVBLGT(KK+ABASE,0,AW4) AW1=AW1-AW2*AW3/AW4 CALL XVBLST(IJ+ABASE,0,AW1) 70 CONTINUE DO 80 J=1,M IJ=(I-1)*RRW+J KJ=(K-1)*RRW+J C80 B(IJ)=B(IJ)-B(KJ)*A(IK)/A(KK) CALL XVBLGT(IJ+BBASE,0,BW1) CALL XVBLGT(KJ+BBASE,0,BW2) BW1=BW1-BW2*AW3/AW4 CALL XVBLST(IJ+BBASE,0,BW1) 80 CONTINUE 90 CONTINUE C SUBSTITUTE BACK NN=(N-1)*RRW+N C NN=N*N CALL XVBLGT(NN+ABASE,0,AW1) IF(AW1.EQ.0.)GOTO 999 DO 110 J=1,M NJ=(N-1)*RRW+J C B(NJ)=B(NJ)/A(NN) CALL XVBLGT(NJ+BBASE,0,BW1) BW1=BW1/AW1 CALL XVBLST(NJ+BBASE,0,BW1) I1MAX=N-1 IF(I1MAX)110,110,95 95 DO 111 I1=1,I1MAX I=N-I1 IJ=(I-1)*RRW+J II=(I-1)*RRW+I I2=I+1 CALL XVBLGT(II+ABASE,0,AW1) DO 100 L=I2,N IL=(I-1)*RRW+L LJ=(L-1)*RRW+J C100 B(IJ)=B(IJ)-A(IL)*B(LJ) CALL XVBLGT(IJ+BBASE,0,BW1) CALL XVBLGT(IL+ABASE,0,AW2) CALL XVBLGT(LJ+BBASE,0,BW2) BW1=BW1-AW2*BW2 CALL XVBLST(IJ+BBASE,0,BW1) 100 CONTINUE C B(IJ)=B(IJ)/A(II) BW1=BW1/AW1 CALL XVBLST(IJ+BBASE,0,BW1) 111 CONTINUE 110 CONTINUE RETURN 999 CONTINUE D=0. RETURN END