CHARACTER*40 INFILE,FILENAME
	LOGICAL EFLG
	real*8 M1(3,4),M2(3,4),M3(3,4),M(3,4),ACC(3,4)
	REAL*8 P1(3),P2(3),T(3)
	character*40 string
	INTEGER*4 MM,NN,I,J

c 	read  operator into m1
c	write(6,*) 'x2= T2 + [rot2] {T1+[rot1]X1}'
5	WRITE(6,*)'ENTER NAME OF rot-trans op (.o):'
	READ(5,1200) INFILE
1200	FORMAT(A)
	IF (INFILE(:5).EQ.'     ') STOP
	  INQUIRE(FILE=INFILE, NAME=FILENAME,EXIST=EFLG)
	  IF (.NOT.EFLG) GOTO 5
        OPEN (UNIT=2,FILE=INFILE,STATUS='OLD',READONLY)
c	READ(2,*) '.SPACE_GROUP_operators    r   12 (3f12.7)'
	READ(2,'(A40)') string
	do 50 j=1,4
C50	read(2,1151)(M1(I,J),I=1,3)
50	read(2,*)(M1(I,J),I=1,3)
1151	FORMAT(3f14.7)
	CLOSE(UNIT=2)
	MM=3
	NN=4
	CALL MATPRN (M1,MM,NN)

	call decomp(M1,M2,NN)

	WRITE(6,*)' RESULTING NCS OPERATOR:'
	WRITE(6,*)' ROTATE MATRIX'
	DO 1100 I=1,3
1100	WRITE(6,9) (M1(I,J),J=1,3)
	WRITE(6,*)' TRANSLATE '
    	WRITE(6,9) (M1(I,4),I=1,3)
9	FORMAT(' ',3F14.6)
	write(6,*)' this operator will be saved as an O datablock RT.O'
	WRITE(6,*)'  Note this operator is for orthogpnal coord, A translation'
	WRITE(6,*)'  You will have to convert for your cell.    '
	WRITE(6,*) 'ENTER NAME FOR output of resulting rot-trans op (.o):'
	READ(5,1200) filename
	open (unit=3,name=FILENAME,status='unknown')
	WRITE(3,*) '.SPACE_GROUP_operators    r   12 (3f14.7)'
	do 1150 j=1,4
1150	write(3,1151)(M1(I,J),I=1,3)
	END

	SUBROUTINE MATPRN(ARRAY,M,N)
	REAL*8 ARRAY(3,4)
	WRITE(6,*)'M, N=',M,N
	DO 10 I=1,M
10	WRITE (6,50) (ARRAY(I,J),J=1,N)
50	format(8f10.5)
  	RETURN
        END

c-----------------------------------------------------------------
	SUBROUTINE DECOMP(R,S1,N)
C	DECOMPOSE GENERAL MATRIX R INTO ORTHONORMAL R AND DIAGONAL S2
	REAL*8 R(3,4),M(3,3),S1(3,4),S2(3,3),S1I(3,3),V(3)
c	test orthonormality of "rotation" matrix:
	write(6,*)
	write(6,*)'Rt*R, test orthonormality:'
	do 510 i=1,3
	do 505 k=1,3
	M(I,K)=0
	do 505 j=1,3
505	M(I,k)=M(I,k)+R(j,i)*R(j,k)		!RT(i,j)*R(j,k)
510	write(6,401) (M(I,k),k=1,3)
401	format(' ',2(3f15.8,3x))

C	FIND S, THE SYMMETRIC? SQRT OF M: SQM*SQM=M, SQM=SQM-INV * M
C	   ITERATIVE METHOD, S1 IS APPROX, S2 WILL BE BETTER APPROX

C	FIRST APPROX: S1 IS DIAGONAL, WITH ELEMENTS = SQRT(DIAG OF M)
	DO 100 I=1,3
	DO  90 J=1,3
90	S1(I,J)=0
100	S1(I,I)=SQRT(M(I,I))

110	CONTINUE

C	COPY S1 TO S1I (= S1-INV) AND INVERT:
	DO 120 I=1,3
	DO 120 J=1,3
120	S1I(I,J)=S1(I,J)
	CALL MATINV(S1I,3,DET)

C	MULTIPLY M BY S1I AND PUT RESULT IN S2
	DO 140 I=1,3
	DO 140 K=1,3
	S2(I,K)=0
	DO 140 J=1,3
140	S2(I,K)=S2(I,K)+S1I(I,J)*M(J,K)

C	MEASURE DIFFERENCE BETWEEN S1 AND S2, AND PUT AVERAGE INTO S2
	SSQE=0
	DO 160 I=1,3
	DO 160 J=1,3
	X=S1(I,J)-S2(I,J)
	SSQE=SSQE+X*X
160	S2(I,J)=S2(I,J)+X/2
	WRITE(6,*)'SSQE=',SSQE

C	IF ERROR OK, GOT TO END AND RETURN. OTHERWISE COPY S2 TO S1 AND REPEAT.
	IF (SSQE.LT.1E-12) GOTO 900
	DO 200 I=1,3
	DO 200 J=1,3
200	S1(I,J)=S2(I,J)
	GOTO 110


C	MULTIPLY R BY S-INV ON THE LEFT, PUT RESULT IN R (temporarily in S2)
900	DO 300 I=1,3
	DO 300 K=1,3
	S2(I,K)=0
	DO 300 J=1,3
300	S2(I,K)=S2(I,K)+R(I,J)*S1I(J,K)
   	
C	TEST IF R*S = ORIG MATRIX (R)
	SSQE=0
	DO 410 I=1,3
	DO 400 K=1,3
	S1I(I,K)=0
	DO 400 J=1,3
400	S1I(I,K)=S1I(I,K)+S2(I,J)*S1(J,K)
	DO 405 K=1,3
	V(K)=R(I,K)-S1I(I,K)
405	SSQE=SSQE+V(K)*V(K)
410	WRITE(6,621)(R(I,K),K=1,3),(S1I(I,K),K=1,3),(V(K),K=1,3)
621	format(' ',4(3f9.4,3x))
	WRITE(6,*)' RMSERROR BETWEEN ORIG MAT AND S*R =',SQRT(SSQE)


C	COPY R-MATRIX FROM S2 TO R. S MATRIX STILL IN S1 (BETTER WAS IN S2, BUT R CALC USING INV OF S1)
	DO 500 I=1,3
	DO 500 J=1,3
500	R(I,J)=S2(I,J)

c	test orthonormality of "rotation" matrix:
	write(6,*)
	write(6,*)'Rt*R, test orthonormality:'
	do 610 i=1,3
	do 605 k=1,3
	M(I,K)=0
	do 605 j=1,3
605	M(I,k)=M(I,k)+R(j,i)*R(j,k)		!RT(i,j)*R(j,k)
610	write(6,401) (M(I,k),k=1,3)

	RETURN
	END

C------------------------------------------------------------------

C------------------------------------------------------------------
	SUBROUTINE MATINV (ARRAY, NORDER, DET)
	real*4 ARRAY, AMAX, SAVE
2	DIMENSION ARRAY(3,3), IK (10), JK (10)
10	DET = 1
11 	DO 100  K=1, NORDER
C
C	FIND LARGEST ELEMENT ARRAY(I,J) IN REST OF MATRIX
C
	AMAX=0
21 	DO 30 I=K, NORDER
     	DO 30 J=K, NORDER
23 	IF (ABS (AMAX) - ABS (ARRAY(I,J) ) ) 24, 24, 30
24 	AMAX=ARRAY(I,J)
     	IK(K) =I
     	JK(K) =J
30 	CONTINUE
C
C INTERCHANGE ROWS AND COLUMNS TO PUT AMAX IN ARRAY(K,K)
C
31 	IF (AMAX) 41, 32, 41
32 	DET=0.
     	GO TO 140
41 	I=IK (K)
     	IF (I-K) 21, 51, 43
43 	DO 50 J=1, NORDER
     	SAVE= ARRAY(K,J)
      	ARRAY(K,J)= ARRAY(I,J)
50 	ARRAY(I,J)= -SAVE
51 	J=JK (K)
     	IF (J-K) 21, 61, 53
53 	DO 60 I=1, NORDER
     	SAVE= ARRAY(I,K)
     	ARRAY(I,K)= ARRAY(I,J)
60 	ARRAY(I,J)= -SAVE
C
C	ACCUMULATE ELEMENTS OF INVERSE MATRIX
C
61 	DO 70 I=1, NORDER
     	IF(I-K) 63, 70, 63
63 	ARRAY(I,K) = -ARRAY(I,K)  / AMAX
70 	CONTINUE
71 	DO 80 I=1, NORDER
     	DO 80 J=1, NORDER
     	IF (I-K) 74, 80, 74
74 	IF (J-K) 75, 80, 75
75 	ARRAY(I,J) = ARRAY(I,J) + ARRAY(I,K) *ARRAY(K,J)
80 	CONTINUE
81 	DO 90 J=1, NORDER
     	IF (J-K) 83, 90, 83
83 	ARRAY(K,J) = ARRAY(K,J) / AMAX
90 	CONTINUE
     	ARRAY(K,K) = 1. / AMAX
c100 	DET=DET *AMAX
100	continue
C
C 	RESTORE ORDERING OF MATRIX
C 
101 	DO 130 L=1, NORDER
       	K= NORDER  - L +1
       	J= IK (K)
       	IF (J-K) 111, 111, 105
105 	DO 110  I=1, NORDER
       	SAVE= ARRAY(I,K)
        ARRAY(I,K)= -ARRAY(I,J)
110 	ARRAY(I,J)= SAVE
111 	I= JK(K)
       	IF (I-K) 130, 130, 113
113 	DO 120  J=1, NORDER
       	SAVE= ARRAY(K,J)
       	ARRAY(K,J) = -ARRAY(I,J)
120 	ARRAY(I,J)= SAVE
130 	CONTINUE
140  	RETURN
	END