GRFMatrixMath

//+------------------------------------------------------------------+
//|                                                   MatrixMath.mq4 |
//|                                    Copyleft© 2007, GammaRatForex |
//|                                       http://gwmmarat.com/Forex/ |
//+------------------------------------------------------------------+
#property copyright "Copyleft© 2007, GammaRatForex"
#property link      "http://gwmmarat.com/Forex/"
#property library

//+------------------------------------------------------------------+
//| MatrixAdd                                                      |
//+------------------------------------------------------------------+
double determinant_fixed(double a[][], int k){
   double z;
    switch(k){
      case 1: return(a[0][0]);
         break;
      case 2: return( a[0][0]*a[1][1]-a[1][0]*a[0][1]);
         break;
      case 3:
         z = a[0][0]*a[1][1]*a[2][2] + a[0][1]*a[1][2]*a[2][0] + 
             a[0][2]*a[1][0]*a[2][1] -
            (a[2][0]*a[1][1]*a[0][2] + a[1][0]*a[0][1]*a[2][2] + 
             a[0][0]*a[2][1]*a[1][2]);
         return(z);
         break;
      default:      
          Print("array_range too large for determinant calculation");
         return(0);
   }  
}
double determinant(double a[][]){
   double z;
   if(ArrayRange(a,0) != ArrayRange(a,1)){
      Print("Non-square array entered into determinant");
      return(0);
    }
    switch(ArrayRange(a,0)){
      case 1: return(a[0][0]);
         break;
      case 2: return( a[0][0]*a[1][1]-a[1][0]*a[0][1]);
         break;
      case 3:
         z = a[0][0]*a[1][1]*a[2][2] + a[0][1]*a[1][2]*a[2][0] + 
             a[0][2]*a[1][0]*a[2][1] -
            (a[2][0]*a[1][1]*a[0][2] + a[1][0]*a[0][1]*a[2][2] + 
             a[0][0]*a[2][1]*a[1][2]);
         return(z);
         break;
      default:      
          Print("array_range too large for determinant calculation");
         return(0);
   }  
}
  
int MatrixAdd(double a[][], double b[][],double & c[][],double w=1){
   int i,j;
   for(i=0;i<ArrayRange(c,0);i++){
      for(j=0;j<ArrayRange(c,1);j++){
         c[i][j] = a[i][j]+w*b[i][j];
      }
   }
    return(0);   
}
int MatrixMul(double a[][], double b[][],double & c[][]){
   int i,j, k;
   if(ArrayRange(a,1) != ArrayRange(b,0)){
      Print("array Range Mismatch in Matrix Mul");
      return(-1);
   }
   for(i=0;i<ArrayRange(a,0);i++){
      for(j=0;j<ArrayRange(b,1);j++){
         c[i][j] = 0;
         for(k=0;k<ArrayRange(a,1);k++){
           c[i][j] += (a[i][k]*b[k][j]);
         }  
      }
   }
  return(0);   
}
int MatrixScalarMul(double b,double & a[][]){
   int i,j;
   for(i =0;i<ArrayRange(a,0);i++){
      for(j=0;j<ArrayRange(a,1);j++){
         a[i][j] *= b;
      }
   }
}         
int MatrixInvert(double a[][],double & b[][]){
   double d;
   double temp1[4][4],temp2[4][4];
   int i,j,i1,j1,k;
   
   if(ArrayRange(a,0) != ArrayRange(a,1)){
      Print("Matrix Inverse attempted on non square matrix. Not implemented yet");
      return(-1);
   }
   if(ArrayRange(a,1) > 3){
      Print("Bugger off, not implemented yet");
      return(-1);
   }
   d = determinant(a);
   if(MathAbs(d)<0.000001){
    Print("unstable matrix");
    d = 1;
   }
   k = ArrayRange(a,0);
   if(k==0){
      b[0][0]=1./a[0][0];
      return(0);
   }
   d = determinant(a);
   for(i1=0;i1<k;i1++){
      for(j1=0;j1<k;j1++){
         temp1[i1][j1] = a[i1][j1];
      }
   }   
   for(i=0;i<k;i++){
      for(j=0;j<k;j++){
         //copy the arry. There must be a better way!!
         ArrayCopy(temp2,temp1);
         for(j1=0;j1<k;j1++){
            if(j1==j){
               temp2[i][j1] = 1;
            }else{
               temp2[i][j1] = 0;
            }
         }   
         b[i][j]=determinant_fixed(temp2,k)/d;
      }
   }
}   
int MatrixTranspose(double a[][],double & b[][]){
   int i,j;
   for(i=0;i<ArrayRange(a,0);i++){
      for(j=0;j<ArrayRange(a,1);j++){
         b[j][i] = a[i][j];
      }
   }
   return(0);
 }
int MatrixZero(double & a[][]){
   int i,j; 
   for(i= 0;i<ArrayRange(a,0);i++){
      for(j=0;j<ArrayRange(a,1);j++){
         a[i][j]=0;
      }
   }
   return(0);
} 
int MatrixEye(double & a[][]){
   int i,j; 
   for(i= 0;i<ArrayRange(a,0);i++){
      for(j=0;j<ArrayRange(a,1);j++){
         a[i][j]=0;
         if(i==j)a[i][j]=1;
      }
   }
   return(0);
} 
int MatrixPrint(double a[][]){
   int i,j; 
   for(i= 0;i<ArrayRange(a,0);i++){
      for(j=0;j<ArrayRange(a,1);j++){
      Print(i," ,",j," ,",a[i][j]);
      }
   }
   return(0);
}                                     
//+------------------------------------------------------------------+