SET NOECHO 
PROC SQRTM
SET NODOECHO NOOUTPUT       
* [AMATRIX]=INPUT MATRIX
* [XMATRIX]=OUTPUT MATRIX
* Procedure for obtaining the square root of a matrix
MATRIX A_=[AMATRIX]
GEN1 N_=$ROWS
MATRIX X_=IDEN(N_)
* Allow up to 20 iterations (but we probably don't need that many)
DO #=1,20
MATRIX X_=.5*(X_+A_*INV(X_))
* See how close we are by first turning the matrix into a vector
MATRIX C_=VEC(X_*X_)-VEC(A_)
GEN1 NC_=$ROWS    
* Get the sum of squared discrepancies
MATRIX CC_=C_'C_
* If we have converged then stop the iterations
SAMPLE 1 NC_
ENDIF(C_.LT.1E-15)
ENDO
ATITLE_: Number of iterations 
PRINT ATITLE_ $DO / NONAME 
MATRIX [XMATRIX]=X_
DELETE / ALL_
SET DOECHO OUTPUT      
PROCEND
SET ECHO