/********************************************************************/
/*                                                                  */
/*  Module		: Matrix calculus                               */
/*                                                                  */
/*  Version 	: 1.8                                                */
/*  Last revision : 01/21/93 18:04:56                               */
/*                                                                  */
/*  Description :                                                   */
/*	 Module for matrix multiplication and matrix exponentiation.   */
/*                                                                  */
/*  Functions supplied :                                            */
/*	 VEC matrix_mul ( VEC mat1, VEC mat2 );                        */
/*	 VEC matrix_exp ( VEC mat, int power );                        */
/*                                                                  */
/********************************************************************/

#include "aglobals.h"
#include "fdecla.h"
# include	"storage.h"

extern int dim, dquad;

VEC matrix2_mul ( mat1, mat2 )
VEC mat1, mat2;
{
	register int i, j, d, dq;
	register VEC p, q, r;
	VEC p_mat;
	
	j = d = dim;
	dq = dquad;
	p_mat = CALLOCATE ( dq );
	p = mat1 + dq;
	r = p_mat + dq;
	while ( j-- ) {
		p -= d;
		r -= d;
		q = mat2 + dq;
		i = d;
		while ( i-- ) {
			q -= d;
			if ( *(p+i) )
				add2_vector ( q, r, d );
		}
	}
	return ( p_mat );
}
 
VEC matrixp_mul ( mat1, mat2 )
VEC mat1, mat2;
{
	register int i, j, d, dq;
	register VEC p, q, r;
	register char val;
	VEC p_mat;
	
	j = d = dim;
	dq = dquad;
	p_mat = CALLOCATE ( dq );
	PUSH_STACK();
	p = mat1 + dq;
	r = p_mat + dq;
	while ( j-- ) {
		p -= d;
		r -= d;
		q = mat2 + dq;
		i = d;
		while ( i-- ) {
			q -= d;
			if ( ( val = *(p+i) ) != 0 )
				ADD_MULT ( val, q, r, d );
		}
	}
	POP_STACK();
	return ( p_mat );
}
 
VEC matrix_exp ( mat, power )
VEC mat;
int power;
{
	register int i = 128;
	VEC v_save, result;
	
	v_save = ALLOCATE ( dquad );
	PUSH_STACK();
	result = mat;
	while ( !(power & i ) ) i >>= 1;
	while ( (i >>= 1) != 0 ) {
		result = MATRIX_MUL ( result, result );
		if ( power & i )
			result = MATRIX_MUL ( result, mat );
	}
	copy_vector ( result, v_save, dquad );
	POP_STACK();
	return ( v_save );
}

void
show_mat ( mat )
VEC mat;
{
	register int i, j;
	for ( i = 0; i < dim; i++ ) {
		for ( j = 0; j < dim; j++ )
			printf ( "%1d", mat[i*dim+j] );
		printf ( "\n" );
	}
}

/* end of module matrix calculus */