#ifndef WMMMAT_H
#define WMMMAT_H

/*
 * WMM - Wave-matching method for mode analysis of dielectric waveguides
 * https://wmm.computational-photonics.eu/
 */

/*
 * wmmmat.h
 * linear algebra routines on vectors, matrices represented as 
 * *double, old FORTRAN format ... 
 */

/* set zero: m x n - matrix a */
void matto0(double *a, int m, int n);

/* matrix-multiplication: c = a*b, m x n matrix a, 
				   n x l matrix b,
			      ->   m x l matrix c */
void matmatmult(double *a, int m, int n, double *b, int l, double *c);

/* matrix-copy: b = a, m x n matrices a, b */
void matcopy(double *a, double *b, int m, int n);

/* additive matrix-copy: b = b + a, m x n matrices a, b */
void addmat(double *a, double *b, int m, int n);

/* matrix-comparison: a == b ? 1:ok, 0:no , m x n - matrices */
int matcompare(double *a, double *b, int m, int n);

/* transpose: b = a^t, m x n matrix a -> n x m matrix b*/
void transpose(double *a, int m, int n, double *b);

/* Cholesky-decomposition of a symmetrical, positive matrix  */
/* 0: okay, 1: matrix is not positive numerically */
int choldc(double *mat, int dim, int testm);

/* backward substitution following Cholesky-decomposition */
/* calculate cmat as a solution of lmat*cmat=rmat,   
   m x m - matrix lmat, left-triangular
   m x n - matrix cmat,
   m x n - matrix rmat */
void bwdsubst(double *lmat, int m, double *cmat, int n, double *rmat);

/* backward substitution following Cholesky-decomposition */
/* calculate vector a as a solution of lmat^t * a = x, 
   m x m - matrix lmat, left-triangular, 
   m - dim vector a,
   m - dim vector x */
void vecbwdsubs(double *lmat, int m, double *a, double *x);

/* calculate a solution xvec of the linear system amat*xvec=bvec */
/* amat: symmetrical, real, positive n*n - matrix */
/* only the upper (right) part of amat must be supplied  */
int solvesympos(double *a, int n, double *x, double *b);

/* normalize vector */
double normalizevec(double *v, int n);

/* scalar product of two vectors */
double scalarprod(double *a, double *b, int n);

/* save matrix to file, MATLAB format */
void savemat(double *mat, int m, int n, const char *nam);

/* test of a n x n matrix for symmetry */
void testmatsym(double *mat, int n, double minel);

/* copy dm x dn block at row l, column c from m x n matrix a to b */ 
void extractbloc(double *a, int m, int n, int l, int c, 
		 double *b, int dm, int dn);

/* copy dm x dn matrix b at row l, column c to m x n matrix a */ 
void insertbloc(double *b, int dm, int dn,
       	        double *a, int m, int n, int l, int c); 

/* enforce symmetry */
void restoresymmetry(double *mat, int dim);

/* ------------------------------------------------------------------------ */

/* small eigenvalue / corresp. eigenvector of a real symmetrc matrix */
/* only the upper (right) part has to be supplied */
/* dim: dimension of matrix smat */
/* mode 0: calculate eigenvalue only 
        1: calculate eigenvalue and eigenvector ev */
/* evidx: index of the required eigenvalue, 0: the lowest, 1: the next larger */
/* pos 0: matrix is not guarateed to be positive 
       1: matrix is guarateed to be positive */

double smallesteigval(double *smat, int dim, int mode, double *ev, int evidx, 
		      int pos);

/* ------------------------------------------------------------------------ */

/* all eigenvalues / corresp. eigenvectors of a real symmetrc matrix */
/* only the upper (right) part has to be supplied */
/* dim: dimension of matrix smat */
/* mode 0: calculate eigenvalue only 
        1: calculate eigenvalue and eigenvector ev */
/* based on: "tred2.c", "tqli.c", "dpythag.c"
   (C) Copr. 1986-92 Numerical Recipes Software 5.)13. */

void eigvalnumrec(double *smat, int dim, int mode, double *eval);

#endif // WMMMAT_H