C----------------------------------------------------------------------- C Library of matrix operations for crystal geometry C----------------------------------------------------------------------- SUBROUTINE MATRIX(A,B,C,D,IWHAT) COMMON /IOUASS/ IOUNIT(12) CHARACTER COUT*132 COMMON /IOUASC/ COUT(20) CHARACTER*6 IWHAT DIMENSION A(3,3),B(3,3),C(3,3),D(3,3),E(3,3),V(3) DATA RA/57.29578/ IF (IWHAT .EQ. 'INVERT') GO TO 100 IF (IWHAT .EQ. 'MATMUL') GO TO 120 IF (IWHAT .EQ. 'MATVEC') GO TO 150 IF (IWHAT .EQ. 'VECMAT') GO TO 180 IF (IWHAT .EQ. 'SCALPR') GO TO 210 IF (IWHAT .EQ. 'LENGTH') GO TO 230 IF (IWHAT .EQ. 'ORTHOG') GO TO 250 IF (IWHAT .EQ. 'DETERM') GO TO 270 IF (IWHAT .EQ. 'MVMULT') GO TO 290 IF (IWHAT .EQ. 'VMMULT') GO TO 320 IF (IWHAT .EQ. 'TRNSPS') GO TO 340 IF (IWHAT .EQ. 'SYMOPR') GO TO 370 IF (IWHAT .EQ. 'VECPRD') GO TO 400 IF (IWHAT .EQ. 'COPRIM') GO TO 410 IF (IWHAT .EQ. 'INTRCH') GO TO 440 IF (IWHAT .EQ. 'SUMVEC') GO TO 460 IF (IWHAT .EQ. 'DIFVEC') GO TO 480 IF (IWHAT .EQ. 'COPVEC') GO TO 500 ITP = IOUNIT(6) WRITE (COUT,10000) IWHAT CALL GWRITE (ITP,' ') STOP C----------------------------------------------------------------------- C Invert 3x3 matrix A, put the result in B C----------------------------------------------------------------------- 100 E(1,1) = A(2,2)*A(3,3) - A(2,3)*A(3,2) E(2,1) = -(A(2,1)*A(3,3) - A(2,3)*A(3,1)) E(3,1) = A(2,1)*A(3,2) - A(2,2)*A(3,1) E(1,2) = -(A(1,2)*A(3,3) - A(1,3)*A(3,2)) E(2,2) = A(1,1)*A(3,3) - A(1,3)*A(3,1) E(3,2) = -(A(1,1)*A(3,2) - A(1,2)*A(3,1)) E(1,3) = A(1,2)*A(2,3) - A(1,3)*A(2,2) E(2,3) = -(A(1,1)*A(2,3) - A(1,3)*A(2,1)) E(3,3) = A(1,1)*A(2,2) - A(1,2)*A(2,1) DMAT = A(1,1)*E(1,1) + A(1,2)*E(2,1) + A(1,3)*E(3,1) DO 115 I=1,3 DO 110 J = 1,3 110 B(I,J) = E(I,J)/DMAT 115 CONTINUE GO TO 520 C----------------------------------------------------------------------- C Multiply 3x3 matrices A and B, store result in C C----------------------------------------------------------------------- 120 DO 135 I = 1,3 DO 132 J = 1,3 E(I,J) = 0.0 DO 130 K = 1,3 130 E(I,J) = E(I,J) + A(I,K)*B(K,J) 132 CONTINUE 135 CONTINUE DO 145 I = 1,3 DO 140 J = 1,3 140 C(I,J) = E(I,J) 145 CONTINUE GO TO 520 C----------------------------------------------------------------------- C Multiply matrix A by vector B, store dir. cosines of result in C C----------------------------------------------------------------------- 150 DO 165 I = 1,3 V(I) = 0. DO 160 J = 1,3 160 V(I) = V(I) + A(I,J)*B(J,1) 165 CONTINUE IF(V(1)**2 + V(2)**2 +V(3)**2 .GT. 0) THEN VMOD = SQRT(V(1)**2 + V(2)**2 + V(3)**2) ELSE VMOD = 1 ENDIF DO 170 I = 1,3 170 C(I,1) = V(I)/VMOD GO TO 520 C----------------------------------------------------------------------- C Multiply vector A by matrix B, store dir. cosines of result in C C----------------------------------------------------------------------- 180 DO 195 I = 1,3 V(I) = 0. DO 190 J = 1,3 190 V(I) = V(I) + B(J,I)*A(J,1) 195 CONTINUE VMOD = SQRT(V(1)**2 + V(2)**2 + V(3)**2) DO 200 I = 1,3 200 C(I,1) = V(I)/VMOD GO TO 520 C----------------------------------------------------------------------- C Scalar product of vectors A and B C----------------------------------------------------------------------- 210 S = 0 DO 220 I = 1,3 220 S = S + A(I,1)*B(I,1) C(1,1) = S GO TO 520 C----------------------------------------------------------------------- C length of vector B when A is the metric matrix C----------------------------------------------------------------------- 230 DO 245 I = 1,3 V(I) = 0. DO 240 J = 1,3 240 V(I) = V(I) + A(I,J)*B(J,1) 245 CONTINUE C(1,1) = SQRT(V(1)**2 + V(2)**2 + V(3)**2) GO TO 520 C----------------------------------------------------------------------- C Get the metric matrix C corresponding to cell edges A & angles B C----------------------------------------------------------------------- 250 COSGAS = (COS(B(1,1)/RA)*COS(B(2,1)/RA) - COS(B(3,1)/RA)) COSGAS = COSGAS/(SIN(B(1,1)/RA)*SIN(B(2,1)/RA)) SINGAS = SQRT(1.0 - COSGAS**2) E(1,1) = A(1,1)*SIN(B(2,1)/RA)*SINGAS E(1,2) = 0 E(1,3) = 0 E(2,1) = -A(1,1)*SIN(B(2,1)/RA)*COSGAS E(2,2) = A(2,1)*SIN(B(1,1)/RA) E(2,3) = 0 E(3,1) = A(1,1)*COS(B(2,1)/RA) E(3,2) = A(2,1)*COS(B(1,1)/RA) E(3,3) = A(3,1) DO 265 I = 1,3 DO 260 J = 1,3 260 C(I,J) = E(I,J) 265 CONTINUE GO TO 520 C----------------------------------------------------------------------- C Calculate the determinant D of the vectors A,B,C C----------------------------------------------------------------------- 270 DET = 0. DO 280 I = 1,3 J = I + 1 IF (J .EQ. 4) J = 1 K = 6 - I - J 280 DET = DET + A(I,1)*(B(J,1)*C(K,1) - B(K,1)*C(J,1)) D(1,1) = DET GO TO 520 C----------------------------------------------------------------------- C Multiply matrix A by vector B, store result in C C----------------------------------------------------------------------- 290 DO 305 I = 1,3 E(I,1) = 0 DO 300 J = 1,3 300 E(I,1) = E(I,1) + A(I,J)*B(J,1) 305 CONTINUE DO 310 I = 1,3 310 C(I,1) = E(I,1) GO TO 520 C----------------------------------------------------------------------- C Multiply vector A by matrix B, store result in C C----------------------------------------------------------------------- 320 DO 335 I = 1,3 C(I,1) = 0. DO 330 J = 1,3 330 C(I,1) = C(I,1) + A(J,1)*B(J,I) 335 CONTINUE GO TO 520 C----------------------------------------------------------------------- Ctranspose matrix A and put it in B C----------------------------------------------------------------------- 340 DO 355 I = 1,3 DO 350 J = 1,3 350 E(I,J) = A(J,I) 355 CONTINUE DO 365 I = 1,3 DO 360 J = 1,3 360 B(I,J) = E(I,J) 365 CONTINUE GO TO 520 C----------------------------------------------------------------------- C Get the symmetry-equivalent of an atom C----------------------------------------------------------------------- 370 DO 390 I = 1,3 C(I,1) = 0. DO 380 J = 1,3 380 C(I,1) = C(I,1) + A(I,J)*B(J,1) J = 4 390 C(I,1) = C(I,1) + A(I,J)/12. GO TO 520 C----------------------------------------------------------------------- C Vector product C = A x B C----------------------------------------------------------------------- 400 C(1,1) = A(2,1)*B(3,1) - A(3,1)*B(2,1) C(2,1) = A(3,1)*B(1,1) - A(1,1)*B(3,1) C(3,1) = A(1,1)*B(2,1) - A(2,1)*B(1,1) GO TO 520 C----------------------------------------------------------------------- C Make coprime integers (the smallest non-zero integer will be 1) C----------------------------------------------------------------------- 410 SMALL = 2. DO 420 I = 1,3 IF (ABS(A(I,1)) .LE. 0.1 .OR. ABS(A(I,1)) .GE. SMALL) GO TO 420 SMALL = ABS(A(I,1)) 420 CONTINUE DO 430 I = 1,3 INDEX = A(I,1)/SMALL + 0.5 IF (A(I,1) .LT. 0.) INDEX = A(I,1)/SMALL - 0.5 430 B(I,1) = INDEX GO TO 520 C----------------------------------------------------------------------- C Interchange two vectors A and B C----------------------------------------------------------------------- 440 DO 450 I = 1,3 SAVE = A(I,1) A(I,1) = B(I,1) 450 B(I,1) = SAVE GO TO 520 C----------------------------------------------------------------------- C Sum of vectors C = A + B C----------------------------------------------------------------------- 460 DO 470 I = 1,3 470 C(I,1) = A(I,1) + B(I,1) GO TO 520 C----------------------------------------------------------------------- C Vector difference C = A - B C----------------------------------------------------------------------- 480 DO 490 I = 1,3 490 C(I,1) = A(I,1) - B(I,1) GO TO 520 C----------------------------------------------------------------------- C Vector copy B = A C----------------------------------------------------------------------- 500 DO 510 I = 1,3 510 B(I,1) = A(I,1) GO TO 520 520 RETURN 10000 FORMAT(' Matrix operation ',A6,' is not programmed') END