	Program MTEST		!Test package for MATRIX

	real*4	A0(40),A(40),B(40),C(40),D(40)
	real*4	X0(10),W(10),X(10),Y(10),Z(10)
	integer*2	i,j,m,n,NF(10)
	data		m,n/3,3/


100	type 1000
1000	format('0','Test program for MATRIX package',/
	1 x,'Please enter m, n (2-6, m>n): ',$)
	accept 1005,m,n
1005	format(2I3)
	if(m.lt.2.or.m.gt.6)goto 100
	if(n.lt.2.or.n.gt.6)goto 100
	if(m*n.gt.40)goto 100

	open(unit=1,name='MATRIX.DAT',type='OLD')
	read(1,2000) A0
	read(1,2000) X0
2000	format(10F8.4)
	close(unit=1)

	type 1010
1010	format('0','MXWRT    A: 1.0 to 12.0')
	call mxwrt(A0,m,n)

	type 1020
1020	format('0','MXMOV    C: same as A')
	call mxmov(A0,m,n,C)
	call mxwrt(C,m,n)

	type 1030
1030	format('0','MXUNT    I')
	call mxunt(B,m,n)
	call mxwrt(B,m,n)

	type 1040
1040	format('0','MXZER    0')
	call mxzer(B,m,n)
	call mxwrt(B,m,n)

	type 1050
1050	format('0','MXSCM    B = 2*A')
	call mxmov(A0,m,n,B)
	call mxscm(B,m,n,2.0)
	call mxwrt(B,m,n)

	type 1060
1060	format('0','MXADD    D = A + C = 2*A')
	call mxadd(A0,C,m,n,D)
	call mxwrt(D,m,n)

	type 1070
1070	format('0','MXSUB    C = D - B = 2*A-2*A = 0')
	call mxsub(D,B,m,n,C)
	call mxwrt(C,m,n)

	type 1080
1080	format('0','MXTRAP   B = A''')
	call mxmov(A0,m,n,B)
	call mxtrap(B,m,n)
	call mxwrt(B,n,m)

	type 1090
1090	format('0','MXNORM')
	call mxnorm(A0,m,n,s)
	type 1100,s
1100	format(x,'Frobenius norm of A is ',F10.4)

	type 1110
1110	format('0','MXMLT    D =  A * (A'')')
	call mxmov(A0,m,n,B)
	call mxtrap(B,m,n)
	call mxmlt(A0,B,m,n,m,D)
	call mxwrt(D,m,m)

	type 1120
1120	format('0','MXMBT    C =  A * A''')
	call mxmov(A0,m,n,B)
	call mxmbt(A0,B,m,n,m,C)
	call mxwrt(C,m,m)

	type 1130
1130	format('0','Check:   B = D - C = 0')
	call mxsub(D,C,m,n,B)
	call mxwrt(B,m,m)

	type 1135
1135	format('0','MXMAT   C = A'' * A')
	call mxmov(A0,m,n,B)
	call mxmat(B,A0,n,m,n,C)
	call mxwrt(C,n,n)

	type 1136
1136	format('0','Check:   (A'')*A - A''*A should = 0')
	call mxtrap(B,m,n)
	call mxmlt(B,A0,n,m,n,D)
	call mxsub(D,C,n,n,B)
	call mxwrt(B,n,n)


	type 1140
1140	format('0','GINV     A.X = C:   find 1/A, & X')
	type 1150
1150	format(x,'Matrix A:')
	call mxmov(A0,m,n,A)
	call mxwrt(A,m,n)
	type 1160
1160	format('0','Vector C:')
	call mxmov(X0,n,1,Y)
	call mxwrt(Y,1,n)

	type *
	call ginv(A,m,n,D,Z,NF,Y,X,2,W)
	type 1170
1170	format('0','Inverse 1/A:')
	call mxwrt(A,n,m)

	type 1175
1175	format('0','Trimmed inverse 1/A:')
	call mxtrim(A,n,m,1.0e-6)
	call mxwrt(A,n,n)

	type 1180
1180	format('0','Solution vector X:')
	call mxwrt(X,n,1)

	type 1185
1185	format('0','Trimmed solution vector X:')
	call mxtrim(X,n,1,1.0e-6)
	call mxwrt(X,n,1)

	type 1190
1190	format('0','Check: B = 1/A * A, which should = I within roundoff')
	call mxmlt(A,A0,n,m,n,B)
	call mxwrt(B,n,n)

	type 1200
1200	format('0','Trimmed B, which should = I exactly')
	call mxtrim(B,n,n,1.0e-5)
	call mxwrt(B,n,n)

	type 1210
1210	format('0','Check equation: A.X, which should = C within roundoff')
	call mxmlt(A0,X,m,n,1,C)
	call mxwrt(C,m,1)

	type 1215
1215	format('0','Trimmed difference A.X - C, which should = 0')
	call mxsub(X0,C,n,1,B)
	call mxtrim(B,n,1,1.0e-5)
	call mxwrt(B,n,1)

	type 1220
1220	format('0','SVD: A = U.Q.V'' ')
	call mxmov(A0,m,n,A)
	type 1230
1230	format('0','Matrix A:')
	call mxwrt(A,m,n)

	type 1240
1240	format('0','SVD itself')
	call SVD(A,m,n,B,C,D)

	type 1250
1250	format('0','Matrix U')
	call mxwrt(B,m,m)

	type 1260
1260	format('0','Check: U.U'' =I')
	call mxmbt(B,B,m,m,m,A)
	call mxtrim(A,m,m,1.0e-5)
	call mxwrt(A,m,m)

	type 1270
1270	format('0','Check: V''.V =I')
	call mxmbt(C,C,n,n,n,A)
	call mxtrim(A,n,n,1.0e-5)
	call mxwrt(A,n,n)

	type 1280
1280	format('0','Matrix Q: Singular values of A''.A, ordered')
	call mxwrt(D,m,n)

	type *
	stop
	end
