## -*-SM-*- Set SM-mode in emacs
## Support for square 2-d matrices
##
## a matrix is defined as a collection of vectors
## as a_0[] a_1[] a_2[] etc (supposed square ...)
## The name of the matrix is given without the `_' , so
## minv a inverts the matrix given by a_0[] a_1[] a_2[] ...
matrix		# Introduction to matrix utilities
		echo a matrix is defined as a collection of vectors
		echo as a_0[] a_1[] a_2[] etc (supposed square ...)
		echo The name of the matrix is given without the `_' , so
		echo minv a inverts the matrix given by a_0[] a_1[] a_2[] ...
		echo
		echo   mprint A                 print A
		echo   mdimen A                 print dimension of A
		echo   mdeclare A N             make A an NxN matrix
		echo   mset A B                 A = B
		echo   madd A B C               A = B+C
		echo   minv A                   A = A^{-1}
		echo   mmult A B C              A = B*C
		echo   mmult_c A B c            A = B*c  (c a scalar)
		echo   mmult_v u A v            u = A*v  (v a vector)
		echo   mmult_vT u v A           u = v^T*A (v a vector)
		echo   mtrans A B               A = B^T
		echo   vprint v                 print v
mdeclare 2	# declare $1 to be a $2x$2 square matrix	
		do 3=0,$2 - 1{
		   set dimen($1_$3) = $2
		}
mdimen	1	# print matrix $1's dimension
		echo $1: $(dimen($1_0))
pmatrix		## alias for mprint
		mprint
mprint	1	# print the matrix $1  (given as $1_0 $1_1 etc -- see minv)
		define 2 (dimen($1_0))
		define 3 < $1_0>
		define 4 <%11.6g>
		if( $2 > 1 ) {
		   do 5=1,($2-1) {
		      define 3 <$3 $1_$5> 
		      define 4 <$4 %11.6g>
		   }
		}
		local define print_noheader 1
		print '$4\n' { $!!3 }
		#
madd	3	# set matrix $1 equal to $2+$3
		mdeclare $1 dimen($2_0)
		do 4=0,dimen($2_0) - 1 {
		   set $1_$4 = $2_$4+$3_$4
		}
qminv  		## matrix inverse, an alias for minv
		minv
minv	1	# Quick matrix inversion, done in place.
		# usage: minv matrix, where
		# $1 is the matrix to be inverted (replaced by the inverse)
		define _n_ (dimen($1_0) - 1)
		do 2=0,$_n_ {
		   set _x_ = $1_$2[$2]
		   set $1_$2[$2] = 1.0
		   set $1_$2 = $1_$2 / _x_
		   do 3=0,$_n_ {
		      if( $2 != $3 ) {
		         set _x_ = $1_$3[$2]
		         set $1_$3[$2] = 0.0
		         set $1_$3 = $1_$3 - _x_ * $1_$2
		      }
		   }
		}
mset	2	# set matrix $1 equal to $2
		do 3=0,dimen($2_0) - 1 {
		   set $1_$3 = $2_$3
		}
mmult	3	# set matrix $1 equal to $2*$3
		local define _n dimen($2_0)
		set _tmp local  set dimen(_tmp) = $_n
		mdeclare $1 $_n
		do 4=0,$_n - 1{
		   do 5=0,$_n - 1{
		      set _tmp[$5] = $2_$5[$4]
		   }
		   do 5=0,$_n - 1{
		      set $1_$5[$4] = SUM(_tmp*$3_$5)
		   }
		}
		#
mmult_c	3	# set matrix $1 equal to $2*$3 where $3 is a scalar
		mdeclare $1 dimen($2_0)
		do 4=0,dimen($1_0) - 1 {
		   set $1_$4 = $2_$4*$3
		}
mmult_v	3	# set vector $1 equal to $2*$3 where $3 is a vector
		define _n dimen($2_0)
		set dimen($1) = $_n
		set dimen(_tmp) = $_n
		do 4=0,$_n - 1{
		   do 5=0,$_n - 1{
		      set _tmp[$5] = $2_$5[$4]
		   }
		   set $1[$4] = SUM($3*_tmp)
		}
mmult_vT 3	# set vector $1 equal to $2^T*$3 where $2 is a vector
		define _n dimen($3_0)
		set dimen($1) = $_n
		do 4=0,$_n - 1{
		   set $1[$4] = SUM($2*$3_$4)
		}
mtrans	2	# set matrix $1 equal to transpose of $2
		mdeclare $1 dimen($2_0)
		define i local  define j local
		do i=0,dimen($2_0) - 1 {
		   do j=0,dimen($2_0) - 1 {
		      set $1_$i[$j] = $2_$j[$i]
		   }
		}
		#
vprint	1	# print the vector $1
		local define print_noheader 1
		print '%11.6g\n' < $1 >
		#