C========================================================== C KEFALAIO 3. SYSTHMATA GRAMMIKWN EJISWSEWN C========================================================== C C PARADEIGMA KYRIWS PROGRAMMATOS GIA TIS DECOMP KAI SOLVE C INTEGER NDIM PARAMETER (NDIM = 10) C REAL A(NDIM,NDIM), B(NDIM), WORK(NDIM), COND, CONDP1 INTEGER IPVT(NDIM), I, N C N = 3 C A(1,1) = 10. A(2,1) = -3. A(3,1) = 5. A(1,2) = -7. A(2,2) = 2. A(3,2) = -1. A(1,3) = 0. A(2,3) = 6. A(3,3) = 5. C CALL DECOMP(NDIM, N, A, COND, IPVT, WORK) WRITE(6,3) COND C CONDP1 = COND + 1 IF (CONDP1 .EQ. COND) THEN WRITE(6,4) STOP ENDIF C B(1) = 7. B(2) = 4. B(3) = 6. C CALL SOLVE(NDIM, N, A, B, IPVT) WRITE(6,5) WRITE(6,6) (B(I), I = 1, N) C 3 FORMAT(/,1X,'Condition number =', E16.9,/) 4 FORMAT(1X,'Matrix is practically singular') 5 FORMAT(1X,'Solution:') 6 FORMAT(1X,10F10.5) C STOP END C C------------------------------------------------------- C SUBROUTINE DECOMP(NDIM, N, A, COND, IPVT, WORK) INTEGER NDIM, N REAL A(NDIM,N), COND, WORK(N) INTEGER IPVT(N) C C ΤΡΙΓΩΝΟΠΟΙΕΙ ΕΝΑΝ ΠΡΑΓΜΑΤΙΚΟ (REAL) ΠΙΝΑΚΑ ΜΕ ΑΠΑΛΟΙΦΗ GAUSS C KAI EKTIMA THN ΚΑΤΑΣΤΑΣΗ ΤΟΥ ΣΥΣΤΗΜΑΤΟΣ. C C ΧΡΗΣΙΜΟΠΟΙΗΣΤΕ ΤΗΝ SOLVE ΓΙΑ ΝΑ ΥΠΟΛΟΓΙΣΕΤΕ ΤΗΝ ΛΥΣΗ ΤΟΥ C ΓΡΑΜΜΙΚΟΥ ΣΥΣΤΗΜΑΤΟΣ. C C ΕΙΣΟΔΟΣ... C C NDIM = ΔΗΛΩΜΕΝΗ ΔΙΑΣΤΑΣΗ ΤΟΥ ARRAY ΠΟΥ ΠΕΡΙΕΧΕΙ ΤΟΝ ΠΙΝΑΚΑ Α. C N = ΤΑΞΗ ΤΟΥ ΠΙΝΑΚΑ. C A = ΠΙΝΑΚΑΣ ΠΟΥ ΘΑ ΤΡΙΓΩΝΟΠΟΙΗΘΕΙ. C C ΕΞΟΔΟΣ... C C A = Ο Α ΠΕΡΙΕΧΕΙ ΕΝΑΝ ΤΡΙΓΩΝΙΚΟ ΠΙΝΑΚΑ U ΚΑΙ ΕΝΑΝ KAΤΩ C ΤΡΙΓΩΝΙΚΟ ΠΙΝΑΚΑ L ΤΕΤΟΙΟΥΣ ΩΣΤΕ: C (ΠΙΝΑΚΑΣ ΑΝΤΙΜΕΤΑΘΕΣΗΣ) * Α = L * U . C C COND = ΕΚΤΙΜΗΣΗ ΤΟΥ ΔΕΙΚΤΗ ΚΑΤΑΣΤΑΣΗΣ ΤΟΥ Α. ΓΙΑ ΤΟ ΓΡΑΜΜΙΚΟ C ΣΥΣΤΗΜΑ A * X = B , ΜΕΤΑΒΟΛΕΣ ΣΤΟΝ Α ΚΑΙ ΣΤΟ Β ΕΙΝΑΙ C ΠΙΘΑΝΟN ΝΑ ΠΡΟΚΑΛΕΣΟΥΝ COND ΦΟΡΕΣ ΜΕΓΑΛΥΤΕΡΕΣ ΔΙΑΤΑ- C ΡΑΧΕΣ ΣΤΟ Χ. C ΑΝ COND + 1.0 = COND O Α ΕΙΝΑΙ ΠΡΑΚΤΙΚΑ ΜΗ ΑΝΤΙΣΤΡΕ- C ΨΙΜΟΣ (ΩΣ ΠΡΟΣ ΤΗΝ ΧΡΗΣΙΜΟΠΟΙΟΥΜΕΝΗ ΑΚΡΙΒΕΙΑ). C Ο COND ΠΑΙΡΝΕΙ ΤΗΝ ΤΙΜΗ 1.0Ε+32 ΑΝ O A ΕΙΝΑΙ ΑΚΡΙΒΩΣ C MH ΑΝΤΙΣΤΡΕΨΙΜΟΣ. C C IPVT = ΤΟ ΔΙΑΝΥΣΜΑ ΟΔΗΓΗΣΗΣ. C IPVT(K) = H ΤΑΞΗ ΤΗΣ ΓΡΑΜΜΗΣ ΠΟΥ ΓΙΝΕΤΑΙ ΟΔΗΓΟΣ ΚΑΤΑ C ΤΟ Κ-ΣΤΟ ΒΗΜΑ ΤΗΣ ΑΠΑΛΟΙΦΗΣ. C IPVT(N) = (-1)**(ΠΛΗΘΟΣ ΕΝΑΛΛΑΓΩΝ ΓΡΑΜΜΩΝ). C C ΧΩΡΟΣ ΕΡΓΑΣΙΑΣ... ΤΟ ΔΙΑΝΥΣΜΑ WORK ΠΡΕΠΕΙ ΝΑ ΔΗΛΩΘΕΙ ΣΤΟ C ΚΥΡΙΩΣ ΠΡΟΓΡΑΜΜΑ ΚΑΙ ΝΑ ΣΥΜΠΕΡΙΛΗΦΘΕΙ ΣΤΗΝ ΕΝΤΟΛΗ C ΚΛΗΣΗΣ. ΚΑΤΑ ΤΗΝ ΕΙΣΟΔΟ ΑΓΝΟΟΥΜΕ ΤΟ ΠΕΡΙΕΧΟΜΕΝΟ ΤΟΥ. C ΚΑΤΑ ΤΗΝ ΕΞΟΔΟ ΟΙ ΣΥΝΙΣΤΩΣΕΣ ΤΟΥ ΣΥΝΗΘΩΣ ΔΕΝ C ΕΝΔΙΑΦΕΡΟΥΝ. C C Η ΟΡΙΖΟΥΣΑ ΤΟΥ Α ΜΠΟΡΕΙ ΝΑ ΥΠΟΛΟΓΙΣΘΕΙ ΚΑΤΑ ΤΗΝ ΕΞΟΔΟ ΩΣ C DET(A) = IPVT(N) * A(1,1) * A(2,2) * ... * A(N,N). C REAL EK, T, ANORM, YNORM, ZNORM INTEGER NM1, I, J, K, KP1, KB, KM1, M C IPVT(N) = 1 IF(N .EQ. 1) GO TO 80 NM1 = N - 1 C C ΥΠΟΛΟΓΙΣΜΟΣ ΤΗΣ 1-ΝΟΡΜΑΣ ΤΟΥ Α C ANORM = 0.0 DO 10 J = 1, N T = 0.0 DO 5 I = 1, N T = T + ABS(A(I,J)) 5 CONTINUE IF(T .GT. ANORM) ANORM = T 10 CONTINUE C C ΑΠΑΛΟΙΦΗ GAUSS ME ΜΕΡΙΚΗ ΟΔΗΓΗΣΗ C DO 35 K = 1, NM1 KP1 = K + 1 C C ΕΥΡΕΣΗ ΟΔΗΓΟΥ C M = K DO 15 I = KP1, N IF (ABS(A(I,K)) .GT. ABS(A(M,K))) M = I 15 CONTINUE IPVT(K) = M IF (M .NE. K) IPVT(N) = -IPVT(N) T = A(M,K) A(M,K) = A(K,K) A(K,K) = T C C ΤΟ ΒΗΜΑ ΠΑΡΑΛΕΙΠΕΤΑΙ ΑΝ Ο ΟΔΗΓΟΣ ΕΙΝΑΙ ΜΗΔΕΝ C IF(T .EQ. 0.0) GO TO 35 C C ΥΠΟΛΟΓΙΣΜΟΣ ΠΟΛΛΑΠΛΑΣΙΑΣΤΩΝ C DO 20 I = KP1, N A(I,K) = - A(I,K) / T 20 CONTINUE C C ΕΝΑΛΛΑΓΗ ΚΑΙ ΑΠΑΛΟΙΦΗ ΚΑΤΑ ΣΤΗΛΕΣ C DO 30 J = KP1, N T = A(M,J) A(M,J) = A(K,J) A(K,J) = T IF (T .EQ. 0.0) GO TO 30 DO 25 I = KP1, N A(I,J) = A(I,J) + A(I,K) * T 25 CONTINUE 30 CONTINUE 35 CONTINUE C C COND = (1-ΝΟΡΜΑ ΤΟΥ Α) * (ΕΚΤΙΜΗΣΗ ΤΗΣ 1-ΝΟΡΜΑΣ ΤΟΥ C ΑΝΤΙΣΤΡΟΦΟΥ TOY A) C Η ΕΚΤΙΜΗΣΗ ΓΙΝΕΤΑΙ ΜΕ ΕΝΑ ΒΗΜΑ ΑΝΤΙΣΤΡΟΦΗΣ ΕΠΑΝΑΛΗΨΗΣ ΓΙΑ ΤΟ C `ΜΙΚΡΟ' ΙΔΙΑΖΟΝ ΔΙΑΝΥΣΜΑ. ΑΥΤΟ ΓΙΝΕΤΑΙ ΜΕ ΕΠΙΛΥΣΗ ΔΥΟ ΣΥΣΤΗ- C ΜΑΤΩΝ (ΑΝΑΣΤΡΟΦΟΣ ΤΟΥ Α) * Υ = Ε ΚΑΙ Α * Ζ = Υ, ΟΠΟΥ ΤΟ E C ΕΙΝΑΙ ΕΝΑ ΔΙΑΝΥΣΜΑ ΜΕ ΣΥΝΙΣΤΩΣΕΣ +1 Η' -1 ΠΟΥ ΕΠΙΛΕΓΕΤΑΙ C ΚΑΤΑΛΛΗΛΑ ΓΙΑ ΝΑ ΔΩΣΕΙ ΜΕΓΑΛΟ Υ. Η ΕΚΤΙΜΗΣΗ ΤΟΥ COND ΙΣΟΥΤΑΙ C ΜΕ ΤΟ ΠΗΛΙΚΟ (1-ΝΟΡΜΑ ΤΟΥ Ζ) / (1-ΝΟΡΜΑ ΤΟΥ Υ). C C ΕΠΙΛΥΣΗ ΤΟΥ (ΑΝΑΣΤΡΟΦΟΣ Α) * Υ = Ε C DO 50 K = 1, N T = 0.0 IF (K .EQ. 1) GO TO 45 KM1 = K-1 DO 40 I = 1, KM1 T = T + A(I,K) * WORK(I) 40 CONTINUE 45 EK = 1.0 IF (T .LT. 0.0) EK = -1.0 IF (A(K,K) .EQ. 0.0) GO TO 90 WORK(K) = -(EK+T) / A(K,K) 50 CONTINUE DO 60 KB = 1, NM1 K = N - KB T = 0.0 KP1 = K + 1 DO 55 I = KP1, N T = T + A(I,K) * WORK(K) 55 CONTINUE WORK(K) = T M = IPVT(K) IF (M .EQ. K) GO TO 60 T = WORK(M) WORK(M) = WORK(K) WORK(K) = T 60 CONTINUE YNORM = 0.0 DO 65 I = 1, N YNORM = YNORM + ABS(WORK(I)) 65 CONTINUE C C ΕΠΙΛΥΣΗ ΤΟΥ Α * Ζ = Υ C CALL SOLVE(NDIM, N, A, WORK, IPVT) ZNORM = 0.0 DO 70 I = 1, N ZNORM = ZNORM + ABS(WORK(I)) 70 CONTINUE C C ΕΚΤΙΜΗΣΗ ΤΟΥ ΔΕΙΚΤΗ ΚΑΤΑΣΤΑΣΗΣ C COND = ANORM * ZNORM / YNORM IF (COND .LT. 1.0) COND = 1.0 RETURN C C 1 x 1 C 80 COND = 1.0 IF (A(1,1) .NE. 0.0) RETURN C C ΠΙΝΑΚΑΣ ΑΚΡΙΒΩΣ ΜΗ ΑΝΤΙΣΤΡΕΨΙΜΟΣ C 90 COND = 1.0E+32 RETURN END C C------------------------------------------------------------- C SUBROUTINE SOLVE(NDIM, N, A, B, IPVT) C INTEGER NDIM, N, IPVT(N) REAL A(NDIM,N), B(N) C C ΕΠΙΛΥΣΗ ΤΟΥ ΓΡΑΜΜΙΚΟΥ ΣΥΣΤΗΜΑΤΟΣ A * X = B. C C NA MHN ΧΡΗΣΙΜΟΠΟΙΕΙΤΑΙ ΑΝ Η DECOMP ΕΧΕΙ ΧΑΡΑΚΤΗΡΙΣΕΙ ΤΟΝ C ΠΙΝΑΚΑ Α ΜΗ ΑΝΤΙΣΤΡΕΨΙΜΟ. C C ΕΙΣΟΔΟΣ... C C NDIM = ΔΗΛΩΜΕΝΗ ΔΙΑΣΤΑΣΗ ΤΟΥ ARRAY ΠΟΥ ΠΕΡΙΕΧΕΙ ΤΟΝ ΠΙΝΑΚΑ Α. C N = ΤΑΞΗ ΤΟΥ ΠΙΝΑΚΑ. C A = ΤΡΙΓΩΝΟΠΟΙΗΜΕΝΟΣ ΠΙΝΑΚΑΣ ΤΟΥ ΠΑΙΡΝΟΥΜΕ ΑΠΟ ΤΗΝ DECOMP. C B = ΔΙΑΝΥΣΜΑ ΔΕΥΤΕΡΟΥ ΜΕΛΟΥΣ. C IPVT = ΔΙΑΝΥΣΜΑ ΟΔΗΓΗΣΗΣ ΠΟΥ ΠΑΙΡΝΟΥΜΕ ΑΠΟ ΤΗΝ DECOMP. C C ΕΞΟΔΟΣ... C C B = ΔΙΑΝΥΣΜΑ ΛΥΣΗΣ Χ. C C INTEGER KB, KM1, NM1, KP1, I, K, M REAL T C C ΠΡΑΞΕΙΣ ΑΠΑΛΟΙΦΗΣ ΣΤΟ Β C IF (N .EQ. 1) GO TO 50 NM1 = N - 1 DO 20 K = 1, NM1 KP1 = K + 1 M = IPVT(K) T = B(M) B(M) = B(K) B(K) = T DO 10 I = KP1, N B(I) = B(I) + A(I,K) * T 10 CONTINUE 20 CONTINUE C C ΟΠΙΣΘΟΔΡΟΜΗΣΗ C DO 40 KB = 1, NM1 KM1 = N - KB K = KM1 + 1 B(K) = B(K) / A(K,K) T = - B(K) DO 30 I = 1, KM1 B(I) = B(I) + A(I,K) * T 30 CONTINUE 40 CONTINUE 50 B(1) = B(1) / A(1,1) RETURN END