// this matrix_h01.h file has the matrix storage allocation and deallocation.
// it also has the complex <--> real and small <--> large isomorphisms.




#define matrix_max(a, b) (((a) >= (b))? (a): (b))
#define matrix_min(a, b) (((a) <= (b))? (a): (b))
#define matrix_abs(a) (((a) >= 0)? (a): (-a))
#define matrix_sign(a) (((a) == 0)? (0): ((a) > 0)? (1): (-1))
#define matrix_inv(a) ((matrix_abs(a) <= mat_pref.test_inv)? (0.0): (1.0 / (a)))


MATRIX_DLL_API
double apowerb(double a, unsigned long int b);


//	 WARNING:  neither validation nor allocation is performed by this  subroutine.
//   this subroutine is called by both matrix_propagate_add & matrix_propatate_multiply.
//		hence, it need not be called explicitly if either of the other two subroutines is called.
//   call this subroutine after calling matrix_allocate.  then, modify any elements, as desired.
MATRIX_DLL_API
void matrix_propagate_default(MAT *array_in, MAT *array_out);


//	 WARNING:  neither validation nor allocation is performed by this subroutine.
//   call this subroutine after calling matrix_alloc.  then, modify any elements, as desired.
MATRIX_DLL_API
void matrix_propagate_add(MAT *array_a, MAT *array_b, MAT *array);


//	 WARNING:  neither validation nor allocation is performed by this subroutine.
//   call this subroutine after calling matrix_allocate.  then, modify any elements, as desired.
MATRIX_DLL_API
void matrix_propagate_multiply(MAT *array_a, MAT *array_b, MAT *array);

	
//  the user must call for each matrix, before the termination of his/her program.
MATRIX_DLL_API
int matrix_dealloc(MAT *array);	// safe to call even if not allocated.


//  the library routines allocate storage, as required -- the user need not allocate anything.
//  the array elements imx and jmx have to be pre-assigned.  all other structure elements will be cleared.
//  the array elements  imp, jmp, sym, cmp. & part have to be post-assigned.
//  the array elements inh may be post-assigned.
MATRIX_DLL_API
int matrix_alloc(MAT *array);	// safe to call even if already allocated.


//  returns sum of: 1=invalid storage allocation, 2=not-square symmetric,
//		4=invalid partitioning, 8=invalid complex.
MATRIX_DLL_API
int matrix_validate(MAT *array);


//  returns sum of:  2=incompatible size, 4=incompatible partitioning, 8=incompatible complex.
MATRIX_DLL_API
int matrix_check_compatibility_add(MAT *array_a, MAT *array_b);


//  returns sum of:  2=incompatible size, 4=incompatible partitioning, 8=incompatible complex.
MATRIX_DLL_API
int matrix_check_compatibility_multiply(MAT *array_a, MAT *array_b);


//  the user must call for each matrix, before the termination of his/her program
MATRIX_DLL_API
int matrix_dealloc_cmplx(matrix_cmplx *array);	// safe to call even if not allocated.


//  the library routines allocate storage, as required -- the user need not allocate anything.
//  the array elements imx and jmx have to be pre-assigned.  all other structure elements will be cleared.
//  the array elements  imp, jmp, sym, cmp. & part have to be post-assigned.
//  the array elements inh may be post-assigned.
MATRIX_DLL_API
int matrix_alloc_cmplx(matrix_cmplx *array);	// safe to call even if already allocated.


//  returns sum of: 1=invalid storage allocation, 2=not-square symmetric.
MATRIX_DLL_API
int matrix_validate_cmplx(matrix_cmplx *array);


//  WARNING:  neither validation nor allocation is performed by this subroutine.
MATRIX_DLL_API
void matrix_offset_copy(MAT *array_a, MAT *array_b, unsigned long imx, unsigned long jmx,
				 unsigned long iaoff, unsigned long jaoff, unsigned long iboff,
				 unsigned long jboff);


MATRIX_DLL_API
int matrix_transpose(MAT *array_a, MAT *array);


MATRIX_DLL_API
int matrix_negative(MAT *array_a, MAT *array);


MATRIX_DLL_API
int matrix_half_sum(MAT *array_a, MAT *array_b, MAT *array);


MATRIX_DLL_API
int matrix_half_difference(MAT *array_a, MAT *array_b, MAT *array, MAT *scratch);


//  afterwards, you may want to modify array->sym or array->cmp.  if you set cmp; you have to clear array->part.
MATRIX_DLL_API
int matrix_part_assemble(MAT *array_a, MAT *array_b, MAT *array_c, MAT *array_d, MAT *array); // (A, B; C, D)


MATRIX_DLL_API
int matrix_part_dismantle(MAT *array, MAT *array_a, MAT *array_b, MAT *array_c, MAT *array_d);


//  split to real.  afterwards, you may want to modify array->sym.     A+iB --> (A, B; -B, A)
MATRIX_DLL_API
int matrix_cmplx2_real(MAT *array_a, MAT *array_b, MAT *array_real, MAT *scratch);


//  cmplx to split.  cmplx --> A+iB
MATRIX_DLL_API
int matrix_cmplx_cmplx2(matrix_cmplx *array_cmplx, MAT *array_a, MAT *array_b);


//  cmplx to real.  cmplx --> (A, B; -B, A).
MATRIX_DLL_API
int matrix_cmplx_real(matrix_cmplx *array_cmplx, MAT *array_real,
					  MAT *scratch_1, MAT *scratch_2, MAT *scratch_3);


//  real to split.  (A, B; -B, A) --> A+iB.  afterwards, you may want to modify array->sym.
MATRIX_DLL_API
int matrix_real_cmplx2(MAT *array_real, MAT *array_a, MAT *array_b,
					   MAT *scratch_1, MAT *scratch_2, MAT *scratch_3, MAT *scratch_4);


//  split to cmplx.  A+iB --> cmplx.  will return 32 if the array_a and array_b are not of the same size.
MATRIX_DLL_API
int matrix_cmplx2_cmplx(MAT *array_a, MAT *array_b, matrix_cmplx *array_cmplx);


MATRIX_DLL_API
int matrix_real_cmplx(MAT *array_real, matrix_cmplx *array_cmplx, MAT *scratch_1, MAT *scratch_2, MAT *scratch_3, 
					  MAT *scratch_4, MAT *scratch_5, MAT *scratch_6);

//	finds the trace of a complex matrix.  requires the real-isomorph as its input.
//  will return an error of 4; if the cmp is zero.
MATRIX_DLL_API
int matrix_trace_complex(MAT *array, cmplx *trace);	// use on any matrix whose cmp is equal to 1.


//  if the original trace is not valid; will set this trace and mark it valid.
MATRIX_DLL_API
int matrix_trace(MAT *array, double *trace);	// use on any matrix whose cmp is equal to 0.


MATRIX_DLL_API
int matrix_small_large(MAT *array_a, MAT *array_b, MAT *array_c, MAT *array, MAT *scratch);


MATRIX_DLL_API
int matrix_large_small(MAT *array, MAT *array_a, MAT *array_b, MAT *array_c, 
					   MAT *scratch_1, MAT *scratch_2, MAT *scratch_3);


//  if the array->cmp is one; then, will yield Hermitian and skew-Hermitian, respectively.
MATRIX_DLL_API
int matrix_symmetric_and_skew(MAT *array, MAT *symmetric, MAT *skew_symmetric, MAT *scratch_1, MAT *scratch_2);



MATRIX_DLL_API
int matrix_scalar_multiply(double x, MAT *array_in, MAT *array);


MATRIX_DLL_API
int matrix_copy(MAT *array_in, MAT *array);


MATRIX_DLL_API
int matrix_add(MAT *array_a, MAT *array_b, MAT *array);


MATRIX_DLL_API
int matrix_subtract(MAT *array_a, MAT *array_b, MAT *array);


MATRIX_DLL_API
int matrix_multiply(MAT *array_a, MAT *array_b, MAT *array);




// copyright (c) 2003,4 by R.I. 'Scibor-Marchocki.  last modified Wednesday 07-th April 2004.

// eof
