Inversa matricei nesingulare

Inversarea matricei nesingulare este subprogramul cel mai căutat căci nenumărate aplicații folosesc inversa matricei, de la ecuația de regresie până la rezolvarea de sisteme liniare.

Am tot căutat subprograme pentru inversarea unei matrice și dintre toate am găsit unul interesant în cartea PROBLEME DE ANALIZĂ NUMERICĂ REZOLVATE CU CALCULATORUL avându-i ca autori pe acad. Gheorghe MARINESCU, Irona RIZZOLI, Ileana POPESCU și Cristina ȘTEFAN, carte publicată îmn Editura Academiei RSR, în anul 1987 și care la paginile 171-173 prezintă metoda iterativă HOTELLING-BODEWING de inversare a unei matrice. Se construiește șirul lui NEWTON de torma X(n+1) = X(n) +(I-X(n)*A)*X(n), unde A, I, X, sunt matrice pătrate cu m linii și m coloane, iar cifrele cuprinse între paranteze acolo la indici arată numărul iterației. X(n+1) înseamnă matricea X la iterația n+1, iar X(n) înseamnă matricea X la iterația n. I este matricea unitate, A este matricea ce trebuie inversată. Subprogramul are ca prim pas evaluarea expresiei X1=X0 +(I+X0*A)*X0. În final matricea X1 va conține matricea inversă, iar X0 este matricea de start care se obține prin convenție dintr-o evaluare de forma: X0 = 1./(maxaij) * I. Variabila maxaij este de fapt elementul maxim dintre elementele matricei A. În subprogramul INVMAT, X0 se va referi la matricea de la iterația precedentă și X1 se va referi la matyricea de la iterația următoare.

Se observă că aici trebuie să fie apelate alte subprograme precum:

- subprogram de copiere matrice;

- subprogram de adunare matrice;

- subprogram de scădere matrice;

- subprogram de calcul normă matrice;

- subprogram de înmulțire matrice.

Față de programul dat în cartea pe care am citat-o am făcut o serie de modificări, căci m-am gândit ca utilizatorul să fie scutit în a da el acolo valoarea de start X0, eroarea, criteriu de stopare a procesului iterativ și matricea unitate.

         SUBROUTINE START(A,N,X0)

         MAXAIJ=A(1,1)

         DO  10 I=1,N

         DO  10 J=1,N

         IF(MAXAIJ.GT.A(I,J)) GO TO 10

         MAXAIJ=A(I,J)

10     CONTINUE

         DO  20 I=1,N

         DO  20 J=1,N

20     X0(I,J)=0.

         ALFA=1./MAXAIJ

         DO  30 I=1,N

30     X0(I,I)=ALFA

         RETURN

         END

Această subrutină construiește matricea de start X0 ca fiind matricea unitate înmulțită cu inversul elementului maxim din matricea A.

Pentru scăderea a două matrice pătrate având n linii și n coloane se foolosește subprogramul:

         SUBROUTINE ADMAT(A,B,C,N)

         MAXAIJ=A(1,1)

         DO  10 I=1,N

         DO  10 J=1,N

         C(I,J)=A(I,J)+B(I,J) 

10     CONTINUE

         RETURN

         END

Pentru scăderea a două matrice pătrate având n linii și n coloane se foolosește subprogramul:

         SUBROUTINE ADMAT(A,B,C,N)

         MAXAIJ=A(1,1)

         DO  10 I=1,N

         DO  10 J=1,N

         C(I,J)=A(I,J)-B(I,J) 

10     CONTINUE

         RETURN

         END

Pentru copierea matricei A în matricea C, ambele având n linii și n coloane se foolosește subprogramul:

         SUBROUTINE SCMAT(A,C,N)

         MAXAIJ=A(1,1)

         DO  10 I=1,N

         DO  10 J=1,N

         C(I,J)=A(I,J)

10     CONTINUE

         RETURN

         END

Generarea unei matrice unitate de ordin n se face cu subprogramul:

         SUBROUTINE UNIMT(N,UNIT)

         DO  20 I=1,N

         DO  10 J=1,N

         UNIT(I,J)=0.

10     CONTINUE

         UNIT(I,I)=1.

20     CONTINUE

         RETURN

         END

Norma ANOR a matricei A este calculată cu subprogramul:

         SUBROUTINE NMMAT(A,N,ANOR)

         ANOR=ABS(A1,1))

         DO  20 I=1,N

         DO  10 J=1,N

         IF(ANOR.GT.ABS(A(I,J)) GO TO 10

         ANOR=ABS(A(I,J))

10     CONTINUE

         RETURN

         END

Pentru înmulțirea a două matrice A și B pătrate având n linii și n coloane, rezultatul fiind matricea C tot cu n linii și n coloane se foolosește subprogramul:

         SUBROUTINE PRMAT(A,C,N)

         MAXAIJ=A(1,1)

         DO  20 I=1,N

         DO  20 J=1,N

         CIJ=0.

         DO  10 K=1,N

         CIJ=CIJ+A(I,K)*B(K,J)

10     CONTINUE

         C(I,J)=CIJ

20     CONTINUE

         RETURN

         END

Subprogramul pentru calculul matricei inverse AINV a matricei A cu n linii și n coloane este:

         SUBROUTINE INVMAT(A,AINV,N,IK)

         ITER=0

         Q=.01

         EPS=.0001

         CALL START(A,N,X0)

10     CALL UNIMT(N,U)

         CALL PRMAT(A,X0,X1,N)

         CALL SCMAT(UNIT,X1,DIF,N)

         ITER=ITER+1

         CALL PRMAT(DIF,X0,X1,N)

         CALL ADMAT(X0,X1,X1,N)

         CALL NMMAT(X0,N,ANOR)

         IF(Q**ITER*ANOR/(1.-Q).LE.EPS) GO TO 20

         CALL CPMAT(X1,X0,N)        

         GO TO 10

20     CONTINUE

         CALL CPMAT(X1,AINV,N)

         RETURN

         END

Față de varianta inițială, aici sunt date valorile care dau finețea rezultatului, fără a da flexibilitate utilizatorului să obțină el precizia dorită, din dorința de a nu avea mulți paramertrii în listă.

(27 decembrie 2017)