// File	matrix.cpp
// Date	5.2.1991
// Author	pf


#define TINY 1e-20

#include <stdlib.h>
#include <stdio.h>
#include <string.h>
#include <math.h>
#include <time.h>
#include <iostream.h>
#include "c:\work\matrix\matrix.h"

//-------------

void matrix::error(char *msg1, char *msg2){
	fprintf(stderr,"matrix Fehler: %s  %s\n",msg1,msg2);
	exit(1);
}

//-------------

matrix::matrix(int mrows, int columns, double initval)
{
  p=new matrep;                            // Die matrep-struct erzeugen
  p->r = mrows;
  p->c = columns;

  p->m = new double *[mrows];              // Speicher fuer die Matrixwerte
										   // allokalisieren
  for (int x=0; x<mrows; x++)
	p->m[x] = new double[columns];
  p->n=1;                                  // Refer-Zaehler auf Eins
  for (int i=0; i<mrows; i++){
	for(int j=0; j<columns; j++)
		mval(i,j) = initval;
  }
}

//-------------

matrix::matrix(int mrows, int columns, double *initvalues )
{
//printf("Zeilen = %d, Spalten %d\n",mrows,columns);
  p=new matrep;                            // Die matrep-struct erzeugen
  p->r = mrows;
  p->c = columns;

  p->m = new double *[mrows];              // Speicher fuer Matrixwerte
  for(int x=0; x<mrows; x++)               // reservieren
	p->m[x] = new double[columns];

  p->n = 1;                                // Referenz-Zaehler auf eins
  int c= 0;
  for(int i=0; i<mrows; i++){
	for(int j=0; j<columns; j++){
		mval(i,j)=initvalues[c++];
	}
  }
}

//-------------

matrix::matrix(char *flag, int dimension)  // erzeugt eine Einheitsmatrix
{
  if (flag[0] !='I')
	error("Erzeugen der Einheitsmatrix: "
		  "Matrix(\"I\",dimension)");
  p = new matrep;
  p->r = dimension;
  p->c = dimension;
  p->m = new double *[dimension];
  for(int x=0; x<dimension; x++)
	p->m[x] = new double [dimension];
  p->n = 1;
  for(int i=0; i<dimension; i++){
	for(int j=0; j<dimension; j++)
		mval(i,j) = (i == j ? 1 : 0);
  }
}

//--------------------------------------
// Fehlermeldung, wenn eine Matrixdatei
// nicht im "Standard"-Format ist
//--------------------------------------
//
static char nonstandard[] =
" Datei ist nicht in Standard-Format. Ein Standard-File muss\n"
"mit den Dimensionen der  Matrix beginen, z.B.:\n"
"\t rows=12 columns=14\n oder abgekuerzt :\n\t r=12 c=14\n"
"Zeilen komen vor Spalten und Namen sind in Kleinbuchstaben\n"
"Kommentare werden mit '#' eingeleitet, Daten folgen auf :::\n";

//--------------------------------------

matrix::matrix(char *initfile)             // Matrix aus Datei lesen
{
  const int BSIZE = 120;
  FILE *from;

  if((from = fopen(initfile,"r"))== NULL )
	error("Matrix initializer file kann nicht geoeffnet werden ",
	       initfile);

  char buf[BSIZE], *cp, *cp2;
  int rfound =0, cfound =0, colonsfound =0;// flags
  p= new matrep;

  // --- Parsing des Datei-headers ---
  while(fgets(buf,BSIZE,from))             // fuer jede Header-Zeile
  {					   // Komentare werden mit ANSI-C
					   // strpbrk() entfernt
    if( (cp = strpbrk(buf,"#")) != NULL)   // Komentare suchen
	*cp='\0';			   // Kommentar-Ende-Zeichen
    if( (cp = strpbrk(buf,"r")) != NULL)
      if( (cp2 = strpbrk(buf,"=")) != NULL)
	if ( (cp = strpbrk(buf,"0123456789")) != NULL){
	  p->r = atoi(cp);
	  rfound++;                        // die Flagge zeigt gefudene
					   // row-Wert an
	}

    if( (cp = strpbrk(buf,"c")) != NULL)
      if( (cp2 = strpbrk(cp,"=")) != NULL)
	if ( (cp = strpbrk(cp2,"0123456789")) != NULL){
	  p->c = atoi(cp);
	  cfound++;                       // die Flagge zeigt gefudene
					   // col-Wert an
	}
    if( strstr(buf,":::") != NULL ) {
      colonsfound++;
      break;                               // ... "while"-Schleife verlassen
    }
  }

    if( !rfound || !cfound || !colonsfound) {
      fprintf(stderr,"%s%s",initfile,nonstandard);
      exit(1);
    }

    p->m = new double *[p->r];
    for(int x=0; x < p->r; x++)
      p->m[x] = new double[p->c];

    p->n = 1;                              // Referenz-Zaehler auf 1
    for(int row=0; row < p->r; row++){
      for(int col=0; col < p->c; col++){
	char nb[20];
	fscanf(from,"%s",nb);             // Blank-begrenzte String suchen
	mval(row,col) = atof(nb);         // in Double-Wert kovertieren
	if(ferror(from))
	  error("Problem mit matrix initializer File",initfile);
      }
    }
}

//-------------

matrix::matrix(matrix& x){
  x.p->n++;                                // Referenzzaehler erhoehen
  p = x.p;                                 // Zeiger auf neue matrep
}

//-------------

matrix::~matrix(){
  if (--p -> n == 0){                      // wird Referenz-Zaehler null?
    for (int x=0; x<rows(); x++)
      delete p->m[x];
    delete p->m;
    delete p;
  }
}

//-------------

matrix matrix::operator=(const matrix& rval)
{                                          // aktuellen Wert loeschen
  if(--p->n == 0){                         // kein anderes Objekt benutzt
    for( int x=0; x<rows(); x++)           // den Wert
      delete p->m[x];                      // Matrixwerte loeschen
    delete p->m;                           // matrep loeschen
    delete p;
  }
					   // neuen Wert zuweisen
  rval.p->n++;                             // eine neue Referenz, Referenz-
					   // Zaehler um Eins erholen, pointer
  p = rval.p;                              // umsetzen
  return *this;
}

//-------------

void matrix::write_standard(char *filename, char *msg){
  FILE *to;
  if((to = fopen(filename,"w")) == NULL)
    error("Matrix output file kann nicht geoeffnet werden",filename);
  fprintf(to,"# %s: Matrix file geschrieben nach \"standard\" format\n",
	  filename);
  time_t clock;
  time(&clock);
					   // ------ Header ------
  fprintf(to,"# %s",asctime(localtime(&clock)));
  fprintf(to,"# %s\n",msg);
  fprintf(to,"rows= %d columns= %d\n",rows(),cols());
  fprintf(to,":::\n");
					   // ------ Matrix ------
  for(int row=0; row<rows(); row++){
    for(int col=0; col<cols(); col++){
      fprintf(to,"%6.6g  ",mval(row,col));
      if(ferror(to))
	error("Problem mit matrix output file",filename);
    }
    fprintf(to,"\n");
  }
}

//-------------

matrix matrix::operator+(const matrix& arg){
	if( rows() != arg.rows() || cols() != arg.cols())
    error("Dimensionen muessen gleich sein !");

  matrix sum(rows(),cols());
  for(int i=0; i<rows(); i++){
    for(int j=0; j<cols(); j++)
      sum.mval(i,j)=mval(i,j)+arg.mval(i,j);
  }
  return sum;
}

//-------------

matrix matrix::operator+(const double arg) {
  matrix sum(rows(),cols());
  for(int i=0; i<rows(); i++){
    for(int j=0; j<cols(); j++)
      sum.mval(i,j)=mval(i,j)+arg;
  }
  return sum;
}

//-------------

matrix matrix::operator-(const matrix& arg){
  if( rows() != arg.rows() || cols() != arg.cols())
    error("Dimensionen muessen gleich sein !");

  matrix sum(rows(),cols());
  for(int i=0; i<rows(); i++){
    for(int j=0; j<cols(); j++)
      sum.mval(i,j)=mval(i,j)-arg.mval(i,j);
  }
  return sum;
}

//-------------

matrix matrix::operator-(const double arg) {
  matrix sum(rows(),cols());
  for(int i=0; i<rows(); i++){
    for(int j=0; j<cols(); j++)
      sum.mval(i,j)=mval(i,j)-arg;
  }
  return sum;
}

//-------------

matrix matrix::operator-() {
  matrix unaryminus(rows(),cols());
  for(int i=0; i<rows(); i++){
    for(int j=0; j<cols(); j++)
      unaryminus.mval(i,j)=-mval(i,j);
  }
  return unaryminus;
}

//-------------

matrix matrix::operator*(const matrix& arg) {
  if( cols() != arg.rows())
	error("Spaltenzahl der 2. Matrix muss gleich "
		  "Zeilenzahl der 1. Matrix sein!");

  matrix result(rows(),arg.cols());

  for(int row=0; row<rows(); row++){
	for(int col=0; col<arg.cols(); col++){
	  double sum=0;
	  for(int i=0; i<cols(); i++){
		sum+=mval(row,i)*arg.mval(i,col);
	  }
	  result.mval(row,col)=sum;
	}
  }
  return result;
}

//-------------

matrix matrix::operator*(const double arg) {
  matrix result(rows(),cols());
	for(int i=0; i<rows(); i++)
		for(int j=0; j<cols(); j++)
			result.mval(i,j)=mval(i,j)*arg;
	return result;
}

//-------------

double &matrix::val(int row, int col) {
	if(row >0 && row < rows() && col >0 && col < cols())
		 return (mval(row,col));
	else
		error("Index uzulaessig");
}

//-------------

matrix matrix::transpose() {
	if(rows() != cols())
		error("Matrix muss quadratisch sein!\n");
	matrix trans(rows(), cols());

	for(int row=0; row<rows(); row++)
		for(int col=0; col<cols(); col++)
			trans.mval(col,row)=mval(row,col);
	return trans;
}


//-------------

double matrix::mmin() {
	double temp;
	if(rows() <= 0 || cols() <= 0)
		error("falsche Matrixgroesse fuer min()");
	double minimum = mval(0,0);

	for(int row=0; row<rows(); row++)
		for(int col=0; col<cols(); col++)
			if((temp=mval(row,col))<minimum)
				minimum=temp;
	return minimum;
}

//-------------

double matrix::mmax() {
	double temp;
	if(rows() <= 0 || cols() <= 0)
		error("falsche Matrixgroesse fuer max()");
	double maximum = mval(0,0);

	for(int row=0; row<rows(); row++)
		for(int col=0; col<cols(); col++)
			if((temp=mval(row,col))>maximum)
				maximum=temp;
	return maximum;
}

//-------------

double matrix::mean() {
	double sum=0;
	for(int row=0; row<rows(); row++)
	  for(int col=0; col<cols(); col++)
	    sum+=fabs(mval(row,col));
	return sum/(row*col);
}

//-------------

double matrix::variance(){
	double s_squared = 0;
	double mn = mean();
	for(int row=0; row<rows(); row++)
		for(int col=0; col<cols(); col++) {
			double temp = mval(row,col) - mn;
			temp *= temp;
			s_squared += temp;
		}
	s_squared /= row*col - 1;          // Anzal Elemente - 1
	return s_squared;
}

//-------------

double matrix::determinant(){
  if(rows() != cols() )
    error("Matrix muss quadratisch sein");
  matrix indx(cols());                     // Einen "Index-Vektor" erzeugen
  matrix B(cols());
  int d;

  matrix decomp = lu_decompose(indx,d);
  double determinant = d;
    for(int i=0; i<cols(); i++)
      determinant *= decomp.mval(i,i);
  return determinant;
}

//-------------

matrix matrix::inverse() {
  if (rows() != cols())
    error("Matrix muss quadratisch sein");
  matrix Y("I",rows()); 		   // Einheitsmatrix erzeugen
  matrix indx(cols());                     // "index vector" erzeugen
  matrix B(cols());                        // siehe Press & Flannery
  int d;

  matrix decomp=lu_decompose(indx,d);      // eine Decomposition durchfuehren
  for(int col=0; col<cols(); col++){
    B.copy_column(Y,col,0);
    decomp.lu_back_subst(indx,B);
    Y.copy_column(B,0,col);
  }
  return Y.transpose();
}

//-------------

void matrix::print(char *msg){
  if(*msg) printf("%s:\n",msg);
  for(int row=0; row<rows(); row++){
    for(int col=0; col<cols(); col++)
      printf("% 10.5f ",mval(row,col));
    printf("\n");
    }
}

/************************************************************
  Die privaten Hilsfunktionen fuer determinant() & inverse()
 ************************************************************/

// kopiere from_col von mm to_col von "this"
void matrix::copy_column(matrix &mm, int from_col, int to_col){
  if( rows() != mm.rows())
    error("Zeilenzahl muss gleich sein fuer copy_column()");
  for(int row=0; row<rows();row++)
    mval(row,to_col) = mm.mval(row,from_col);

}

//-------------

void matrix::switch_columns(int col1, int col2){
  matrix temp(rows());
  for(int row=0; row<rows(); row++)
    temp.mval(row,0) = mval(row,col1);     // erste Spalte zwischenspeichern
  for(row=0; row<rows(); row++)
    mval(row,col1) = mval(row,col2);       // col2 nach col1
  for(row=0; row<rows(); row++)
    mval(row,col2) = temp.mval(row,0);     // temp nach col2
}

//-------------

// Eine Echtkopie (fuer die L-U-Dekomposition )
void matrix::deepcopy(matrix & from, matrix & to){
  if(from.rows() != to.rows() || from.cols() != to.cols())
    error("Matrizen muessen gleiche Dimensionen haben fuer deepcopy()");
  for(int row=0; row < from.rows(); row++)
    for(int col=0; col < from.cols(); col++)
      to.mval(row,col) = from.mval(row,col);
}

//-------------

// Matrix skalieren ( fuer L-U-Dekomposition )
matrix matrix::scale() {
  double temp;
  if(rows() <= 0 || cols() <=0)
    error("falsche Matrixgroesse fuer scale()");
  if(rows() != cols())
    error("Matrix muss quadratisch sein fuer scale()");
  matrix scale_vector(rows());
  for(int col=0; col < cols(); col++){
    double maximum =0;
    for(int row = 0; row <rows(); row++)
      if((temp = (double)fabs(mval(row,col))) > maximum)
	maximum=temp;                      // Spalte mit Maximalwert
					   // in dieser Zeile
      if(maximum == 0)
	error("singular matrix in scale()");
      scale_vector.mval(col,0)=1/maximum;  // Skalierungswert
    }
  return scale_vector;
}

//-------------

matrix matrix::lu_decompose(matrix& indx, int& d){
/*
 Liefert die L-U-Dekomposition einer Matrix ( bringt Matrix in Zeilen-
 -Stufen-Form). indx ist ein Output-Vektor. Er beschreibt die Zeilen-
 -Permutationen. d ist ebenfals ein Output-Wert und zeigt an, ob die Anzahl
 der Zielen-Vertauschungen gerade (+1) oder ungerade (-1) waren.
 lu-decompose wird zusammen mit lu-back_subs benutzt, um
 lineare Gleichungssysteme zu ioesen oder die Inverse einer Matrix
 zu errechnen.
*/
  if(rows() != cols())
    error("Matrix muss quadratisch sein fuer L-U decompose!");
  d = 1;                                   // Parity - Check
  int row,col,k,col_max;                   // versch. Zaehler
  double dum;
  double sum;
  double maximum;

  matrix lu_decomp(rows(),cols());    // eine Echtkopie der orig. Matrix
  deepcopy(*this,lu_decomp);
  matrix scale_vector = lu_decomp.scale();

// Crout's Methode:
  for(row = 0; row <rows(); row++){
    if(row >0){
      // Glch. 2.3.12, ausser fuer row=col:
      for(col=0; col<=row-1;col++){
	sum=lu_decomp.mval(row,col);
	if(col > 0){
	  for(k=0;k<=col-1;k++)
	    sum-=lu_decomp.mval(row,k)*lu_decomp.mval(k,col);
	  lu_decomp.mval(row,col) = sum;
	}
      }
    }
  // Initialisierung fuer die Suche nach dem groessten Pivot:
  maximum = 0;
  // i=j siehe Glch. 2.3.12 & i=j+1..N siehe Glch. 2.3.13:
  for(col=row; col<=cols()-1; col++){
    sum=lu_decomp.mval(row,col);
    if(row>0){
      for(k=0; k<=row-1; k++)
	sum-=lu_decomp.mval(k,col) * lu_decomp.mval(row,k);
      lu_decomp.mval(row,col)=sum;
    }
    dum = scale_vector.mval(col,0) * fabs(sum);
    if(dum >= maximum){
      col_max = col;
      maximum = dum;
    }
  }
  // Muessen Zeilen vertauscht werden ?
  if(row != col_max){
    lu_decomp.switch_columns(col_max,row); // Ja ...
    d *=-1;                                // Parity neu setzen
    // Skalierung faktor ebenfals vertauschen :
    dum = scale_vector.mval(col_max,0);
    scale_vector.mval(col_max,0) = scale_vector.mval(col,0);
    scale_vector.mval(row,0) = dum;
  }
  indx.mval(row,0) = col_max;
  // Dividiere zum Schluss das Pivot-Element:
  if(row != rows()-1){
    if(lu_decomp.mval(row,row) == 0.0)
    lu_decomp.mval(row,row) = TINY;
      // falls das Pivot-Element Null ist, ist die
      // Matrix singulaer (entsprechend der Rechengenauigkeit).
      // Fuer manche Anwendungen singulaerer Matrizen ist es
      // einfacher, wenn Null durch TINY ersetzt wird (einen
      // sehr kleinen Wert <> 0 ).
    dum = 1/lu_decomp.mval(row,row);
    for(col=row+1; col <=cols()-1; col++)
      lu_decomp.mval(row,col) *= dum;
  }
 }
  if(lu_decomp.mval(rows()-1,cols()-1) == 0)
    lu_decomp.mval(rows()-1,cols()-1) = TINY;
  return lu_decomp;
}

//-------------

void matrix:: lu_back_subst(matrix & indx, matrix & b) {
/*
  Loest ein System von N linearen Gleichungen A*X=B. "this"
  besteht in diesem Fall aus der LU-dekomponierten Matrix A, die
  durch lu_decompose() erzeugt wird. Indx stammt ebenfalls aus
  lu_decompose() (s.o). B ist der rechtsseitige Vektor der
  Gleichung und enthaelt nach Ende der Funktion die Losung X. Der
  hier verwandte Algorithmus beachtet die Moeglichkeit, dass B mit
  vielen Null-Elementen beginnt, und ist deshalb fuer die Inversion
  einer Matrix geeignet. Siehe Seiten 36-37 in Press & Flannery
*/
  if(rows() != cols())
    error("nichtquadratische lu_decomp Matrix in lu_back_subst()");
  if(rows() != b.rows())
    error("B Vektor falsche Groesse in lu_back_subst()");
  if(rows() != indx.rows())
    error("Indexvektor falsche Groesse in lu_back_subst()");
  int row, col, ll;
  int ii = 0;
  double sum;
  for(col=0; col < cols(); col++){
    ll = (int) indx.mval(col,0);
    sum = b.mval(ll,0);
    b.mval(ll,0) = b.mval(col,0);
    if(ii>=0)
      for(row=ii; row<=col-1; row++)
	sum-=mval(row,col)*b.mval(row,0);
    else
      if(sum !=0)
	ii = col;
    b.mval(col,0) = sum;
  }
  for(col = cols()-1; col >= 0; col--){
    sum = b.mval(col,0);
    if(col < cols()-1)
      for(row = col+1; row <=rows()-1; row++)
	sum -= mval(row,col)*b.mval(row,0);
    //eine Komponente des Losungsvert X speichern
    b.mval(col,0) = sum/mval(col,col);
  }
}

//-----------------------------------------