//===============================================================================
//
//  VectMatr.h
//
//  Library for working with vectors and matrices
//
//  Vertion 0.2
//
//  Designer and developer: Henry Guennadi Levkin
//
//===============================================================================

#include <stdio.h>
#include <math.h>

//-------------------------------------------------------------------------------
#define VECTOR_MAX_DIM  16000000
#define MATRIX_MAX_DIM      4000

//-------------------------------------------------------------------------------
typedef struct
{
  int nD;     // dimension of vector

  double* pV; // array with vector

} SVect;

//-------------------------------------------------------------------------------
typedef struct
{
  int     nD; // dimension of square matrix

  double* pM; // array with matrix (row after row)

  double** m; // array with pointers to rows

} SMatr;

//-------------------------------------------------------------------------------
int VectorConstruct(SVect* pVect, int nDimension)
{
  if(nDimension <= VECTOR_MAX_DIM)
  {
    pVect->nD = nDimension;
  }
  else
  {
    return 1;
  }
  
  pVect->pV = (double*) malloc(nDimension * sizeof(double) );

  if(pVect->pV == NULL) return 2;

  return 0;
}

//-------------------------------------------------------------------------------
int VectorConstructAndFill(SVect* pVect, int nDimension, double* pArr)
{
  int i;

  if(nDimension <= VECTOR_MAX_DIM) 
  {
    pVect->nD = nDimension;
  }
  else
  {
    return 1;
  }

  pVect->pV = (double*) malloc(nDimension * sizeof(double) );
  
  if(pVect->pV == NULL) return 2;
  
  for(i=0; i<nDimension; i++)
  {
    pVect->pV[i] = pArr[i];
  }

  return 0;
}

//-------------------------------------------------------------------------------
int VectorCopyConstruct(SVect* pSrc, SVect* pNew)
{
  int nDim;
  int i; 
  
  nDim = pSrc->nD;

  pNew->pV = (double*) malloc(nDim * sizeof(double) );
  
  if(pNew->pV == NULL) return 1;
  
  pNew->nD = nDim;
  
  for(i=0; i<nDim; i++)
  {
    pNew->pV[i] = pSrc->pV[i];
  }
  
  return 0;
}

//-------------------------------------------------------------------------------
// destructor
void VectorFree(SVect* pVec)
{
  free((char*)(pVec->pV));
}

//-------------------------------------------------------------------------------
int VectorFill(SVect* pVect, double* pArr)
{
  int nDim = pVect->nD;

  memcpy((char*)pVect->pV, (char*)pArr, nDim * sizeof(double) );

  return 0;
}

//-------------------------------------------------------------------------------
int VectorCopy(SVect* pSrc, SVect* pDst)
{
  int i;
  int nDim = pSrc->nD;

  if(pSrc->nD != pDst->nD) return 1;

  for(i=0; i<nDim; i++)
  {
    pDst->pV[i] = pSrc->pV[i];
  }
  
  return 0;
}

//-------------------------------------------------------------------------------
int VectorAddConstant(SVect* pSrc, SVect* pDst, double value)
{
  int i;
  int nDim = pSrc->nD;

  if(pSrc->nD != pDst->nD) return 1;

  for(i=0; i<nDim; i++)
  {
    pDst->pV[i] = pSrc->pV[i] + value;
  }
  
  return 0;
}

//-------------------------------------------------------------------------------
int VectorMultConstant(SVect* pSrc, SVect* pDst, double value)
{
  int i;
  int nDim = pSrc->nD;

  if(pSrc->nD != pDst->nD) return 1;

  for(i=0; i<nDim; i++)
  {
    pDst->pV[i] = pSrc->pV[i] * value;
  }
  
  return 0;
}

//-------------------------------------------------------------------------------
int VectorScalarMult(SVect* pVec1, SVect* pVec2, double* pfResult)
{
  int i;
  double fRes = 0.0;
  int nDim = pVec1->nD;

  if(pVec1->nD != pVec2->nD) return 1;

  for(i=0; i<nDim; i++)
  {
    fRes += ( pVec1->pV[i] * pVec2->pV[i]);
  }
  
  *pfResult = fRes;
  
  return 0;
}

//-------------------------------------------------------------------------------
int VectorAdd(SVect* pVec1, SVect* pVec2, SVect* pResult)
{
  int i;
  int nDim = pVec1->nD;

  if(nDim != pVec2->nD)   return 1;
  if(nDim != pResult->nD) return 2;

  for(i=0; i<nDim; i++)
  {
    pResult->pV[i] = ( pVec1->pV[i] + pVec2->pV[i]);
  }
  
  return 0;
}

//-------------------------------------------------------------------------------
int VectorSubtract(SVect* pVec1, SVect* pVec2, SVect* pResult)
{
  int i;
  int nDim = pVec1->nD;

  if(nDim != pVec2->nD)   return 1;
  if(nDim != pResult->nD) return 2;

  for(i=0; i<nDim; i++)
  {
    pResult->pV[i] = ( pVec1->pV[i] - pVec2->pV[i]);
  }
  
  return 0;
}

//-------------------------------------------------------------------------------
int VectorLength(SVect* pVec, double* pfLength)
{
  int i;
  int nDim = pVec->nD;
  
  double fLength = 0.0;

  for(i=0; i<nDim; i++)
  {
    fLength += ( pVec->pV[i] * pVec->pV[i]);
  }

  *pfLength = sqrt(fLength);

  return 0;
}

//-------------------------------------------------------------------------------
int VectorAbsNorm(SVect* pVec, double* pfNorm)
{
  int i;
  int nDim = pVec->nD;
  
  double fNorm = 0.0;

  for(i=0; i<nDim; i++)
  {
    fNorm += fabs( pVec->pV[i] );
  }

  *pfNorm = fNorm;

  return 0;
}

//-------------------------------------------------------------------------------
int VectorEuclidDist(SVect* pVec1, SVect* pVec2, double* pfDistance)
{
  int i;
  int nDim = pVec1->nD;
  double fDistance = 0.0;
  double fDelta;

  if(pVec1->nD != pVec2->nD) return 1;

  for(i=0; i<nDim; i++)
  {
    fDelta = ( pVec1->pV[i] - pVec2->pV[i]);
    fDistance += (fDelta * fDelta);
  }
  
  *pfDistance = sqrt( fDistance );
  
  return 0;
}

//-------------------------------------------------------------------------------
int VectorChebyshevDist(SVect* pVec1, SVect* pVec2, double* pfDistance)
{
  int i;
  int nDim = pVec1->nD;
  double fDistance = 0.0;
  double fDelta;

  if(pVec1->nD != pVec2->nD) return 1;

  for(i=0; i<nDim; i++)
  {
    fDelta = ( pVec1->pV[i] - pVec2->pV[i]);
    fDistance += fabs( fDelta );
  }
  
  *pfDistance = fDistance;

  return 0;
}

//-------------------------------------------------------------------------------
// angle in redians
int VectorAngle(SVect* pVec1, SVect* pVec2, double* pfAngle)
{
  int i;
  int nDim = pVec1->nD;
  double fScalar = 0.0;
  double fLeng1  = 0.0;
  double fLeng2  = 0.0;
  double fAngle;

  if(pVec1->nD != pVec2->nD) return 1;

  for(i=0; i<nDim; i++)
  {
    fLeng1  += ( pVec1->pV[i] * pVec1->pV[i] );
    fLeng2  += ( pVec2->pV[i] * pVec2->pV[i] );
    fScalar += ( pVec1->pV[i] * pVec2->pV[i] );
  }
  fLeng1 = sqrt(fLeng1);
  fLeng2 = sqrt(fLeng2);
  fAngle = fScalar / fLeng1/ fLeng2;
  
  *pfAngle = acos(fAngle);
  
  return 0;
}

//-------------------------------------------------------------------------------
int VectorPrint(SVect* pVec)
{
  int i;
  int nDim = pVec->nD;

  printf("------------------------------------------------------------------\n");
  printf("Vector dimension = %d\n", nDim);

  for(i=0; i<nDim; i++)
  {
    printf("%.3lf ", pVec->pV[i]);
  }

  printf("\n");

  return 0;
}


//-------------------------------------------------------------------------------
//-------------------------------------------------------------------------------
int MatrixConstruct(SMatr* pMat, int nDim)
{
  int i;

  if(nDim <= MATRIX_MAX_DIM)
  {
    pMat->nD = nDim;
  }
  else
  {
    return 1;
  }
  
  pMat->pM = (double*) malloc(nDim * nDim * sizeof(double) );
  
  if(pMat->pM == NULL) return 2;

  pMat->m = (double**) malloc(nDim * sizeof(double*) );
  
  memset((char*)pMat->pM, 0, nDim * nDim * sizeof(double) );
  
  for(i=0; i<nDim; i++)
  {
    pMat->m[i] = pMat->pM + i * nDim;
  }

  return 0;
}

//-------------------------------------------------------------------------------
int MatrixConstructUnit(SMatr* pMat, int nDim)
{
  int i;

  if(nDim <= MATRIX_MAX_DIM)
  {
    pMat->nD = nDim;
  }
  else
  {
    return 1;
  }
  
  pMat->pM = (double*) malloc(nDim * nDim * sizeof(double) );
  
  if(pMat->pM == NULL) return 2;

  pMat->m = (double**) malloc(nDim * sizeof(double*) );
  
  memset((char*)pMat->pM, 0, nDim * nDim * sizeof(double) );
  
  for(i=0; i<nDim; i++)
  {
    pMat->m[i] = pMat->pM + i * nDim;
  }

  for(i=0; i<nDim; i++)
  {
    pMat->m[i][i] = 1.0;
  }

  return 0;
}
//-------------------------------------------------------------------------------
int MatrixConstructAndFill(SMatr* pMat, int nDim, double* pArr)
{
  int i;
  int nn = nDim * nDim;

  if(nDim <= MATRIX_MAX_DIM)
  {
    pMat->nD = nDim;
  }
  else
  {
    return 1;
  }
  
  pMat->pM = (double*) malloc(nDim * nDim * sizeof(double) );
  
  if(pMat->pM == NULL) return 2;

  pMat->m = (double**) malloc(nDim * sizeof(double*) );
  
  memset((char*)pMat->pM, 0, nDim * nDim * sizeof(double) );
  
  for(i=0; i<nDim; i++)
  {
    pMat->m[i] = pMat->pM + i * nDim;
  }

  for(i=0; i<nn; i++)
  {
    pMat->pM[i] = pArr[i];
  }

  return 0;
}

//-------------------------------------------------------------------------------
int MatrixConstructFromFile(SMatr* pMat, char* pFilename)
{
  int i;
  int nDim;
  int nn;

  FILE* pFile = fopen(pFilename,"r");

  if(pFile==NULL)
  {
    return 1;
  }

  fscanf(pFile,"%d",&nDim);

  if(nDim <= MATRIX_MAX_DIM)
  {
    pMat->nD = nDim;
  }
  else
  {
    return 2;
  }
  
  pMat->pM = (double*) malloc(nDim * nDim * sizeof(double) );
  
  if(pMat->pM == NULL) return 3;

  pMat->m = (double**) malloc(nDim * sizeof(double*) );
  
  for(i=0; i<nDim; i++)
  {
    pMat->m[i] = pMat->pM + i * nDim;
  }
  
  nn = nDim * nDim;

  for(i=0; i<nn; i++)
  {
    fscanf(pFile,"%lf",&(pMat->pM[i]));
  }

  return 0;
}

//-------------------------------------------------------------------------------
int MatrixCopyConstruct(SMatr* pSrc, SMatr* pNew)
{
  int i;
  int nDim = pSrc->nD;
  int nn = nDim * nDim;

  if(nDim <= MATRIX_MAX_DIM)
  {
    pNew->nD = nDim;
  }
  else
  {
    return 1;
  }
  
  pNew->pM = (double*) malloc(nDim * nDim * sizeof(double) );

  if(pNew->pM == NULL) return 2;

  pNew->m = (double**) malloc(nDim * sizeof(double*) );
  
  for(i=0; i<nDim; i++)
  {
    pNew->m[i] = pNew->pM + i * nDim;
  }

  memcpy((char*)pNew->pM, (char*)pSrc->pM, nDim * nDim * sizeof(double) );
  
  return 0;
}

//-------------------------------------------------------------------------------
int MatrixCopy(SMatr* pSrc, SMatr* pDst)
{
  int nDim = pSrc->nD;

  if(pSrc->nD != pDst->nD)
  {
    return 1;
  }

  memcpy((char*)pDst->pM, (char*)pSrc->pM, nDim * nDim * sizeof(double) );

  return 0;
}

//-------------------------------------------------------------------------------
int MatrixAddConstant(SMatr* pSrc, SMatr* pDst, double value)
{
  int i;

  int nn = pSrc->nD * pSrc->nD;

  if(pSrc->nD != pDst->nD)
  {
    return 1;
  }
  
  for(i=0; i<nn; i++)
  {
    pDst->pM[i] = pSrc->pM[i] + value;
  }

  return 0;
}

//-------------------------------------------------------------------------------
int MatrixMultConstant(SMatr* pSrc, SMatr* pDst, double value)
{
  int i;

  int nn = pSrc->nD * pSrc->nD;

  if(pSrc->nD != pDst->nD)
  {
    return 1;
  }
  
  for(i=0; i<nn; i++)
  {
    pDst->pM[i] = pSrc->pM[i] * value;
  }

  return 0;
}

//-------------------------------------------------------------------------------
int MatrixAddMatrix(SMatr* pMat1,SMatr* pMat2, SMatr* pRes)
{
  int i;

  int nn = pMat1->nD * pMat1->nD;

  if(pMat1->nD != pMat2->nD)
  {
    return 1;
  }
  
  for(i=0; i<nn; i++)
  {
    pRes->pM[i] = pMat1->pM[i] + pMat2->pM[i];
  }

  return 0;
}

//-------------------------------------------------------------------------------
int MatrixSubtractMatrix(SMatr* pMat1,SMatr* pMat2, SMatr* pRes)
{
  int i;

  int nn = pMat1->nD * pMat1->nD;

  if(pMat1->nD != pMat2->nD)
  {
    return 1;
  }
  
  for(i=0; i<nn; i++)
  {
    pRes->pM[i] = pMat1->pM[i] + pMat2->pM[i];
  }

  return 0;
}

//-------------------------------------------------------------------------------
int MatrixTranspose(SMatr* pSrc, SMatr* pDst)
{
  int nDim = pSrc->nD;
  int i, j;
  
  if(pSrc->nD != pDst->nD)
  {
    return 1;
  }
  
  for(i=0; i<nDim; i++)
  {
    for(j=0; j<nDim; j++)
    {
      pDst->m[j][i] = pSrc->m[i][j];
    }
  }

  return 0;
}

//-------------------------------------------------------------------------------
int MatrixMultMatrix(SMatr* pMat1,SMatr* pMat2, SMatr* pRes)
{
  int nDim = pMat1->nD;
  int i, j, k;
  
  if(pMat1->nD != pMat2->nD)
  {
    return 1;
  }

  if(pMat1->nD != pRes->nD)
  {
    return 2;
  }

  for(i=0; i<nDim; i++)
  {
    for(j=0; j<nDim; j++)
    {
      pRes->m[i][j] = 0.0;
      for(k=0; k<nDim; k++)
      {
        pRes->m[i][j] += (pMat1->m[i][k] * pMat2->m[k][j]);
      }
    }
  }

  return 0;
}

//-------------------------------------------------------------------------------
int MatrixMultVector(SMatr* pMat, SVect* pVec, SVect* pRes)
{
  int nDim = pMat->nD;
  int i, j;
  
  if(pMat->nD != pVec->nD)
  {
    return 1;
  }

  if(pMat->nD != pRes->nD)
  {
    return 2;
  }

  for(i=0; i<nDim; i++)
  {
    pRes->pV[i] = 0.0;
    for(j=0; j<nDim; j++)
    {
      pRes->pV[i] += (pMat->m[i][j] * pVec->pV[j]);
    }
  }

  return 0;
}

//-------------------------------------------------------------------------------
int MatrixRawsSwap(SMatr* pMat, int nRaw1, int nRaw2, SMatr* pRes)
{
  int nDim = pMat->nD;

  if(pMat->nD != pRes->nD)
  {
    return 1;
  }

  memcpy((char*)pRes->pM, (char*)pMat->pM, nDim * nDim * sizeof(double) );

  memcpy((char*)(pRes->pM + nRaw2*nDim), (char*)(pMat->pM +nRaw1*nDim), nDim * sizeof(double) );

  memcpy((char*)(pRes->pM + nRaw1*nDim), (char*)(pMat->pM +nRaw2*nDim), nDim * sizeof(double) );

  return 0;
}

//-------------------------------------------------------------------------------
int MatrixColumnsSwap(SMatr* pMat, int nColumn1, int nColumn2, SMatr* pRes)
{
  int nDim = pMat->nD;
  int i;
  
  if(pMat->nD != pRes->nD)
  {
    return 1;
  }

  memcpy((char*)pRes->pM, (char*)pMat->pM, nDim * nDim * sizeof(double) );
  
  for(i=0; i<nDim; i++)
  {
    pRes->m[i][nColumn1] = pMat->m[i][nColumn2];
    pRes->m[i][nColumn2] = pMat->m[i][nColumn1];
  }

  return 0;
}

//-------------------------------------------------------------------------------
int MatrixPrint(SMatr* pMat)
{
  int nDim = pMat->nD;
  int i, j;

  printf("===================================================================\n");
  printf("Matrix Dimension = %d\n", nDim);

  for(i=0; i<nDim; i++)
  {
    for(j=0; j<nDim; j++)
    {
      printf("%.4lf ", pMat->m[i][j]);
    }
    printf("\n");
  }
  return 0;
}

//-------------------------------------------------------------------------------
int MatrixPrintToFile(SMatr* pMat, char* pFilename)
{
  int nDim = pMat->nD;
  int i, j;
  
  FILE* pFile = fopen(pFilename,"w");

  fprintf(pFile, "%d\n", nDim);

  for(i=0; i<nDim; i++)
  {
    for(j=0; j<nDim; j++)
    {
      fprintf(pFile, "%.4lf ", pMat->m[i][j]);
    }
    fprintf(pFile, "\n");
  }
  
  fclose(pFile);
 
  return 0;
}

//-------------------------------------------------------------------------------
int MatrixReadFromFile(SMatr* pMat, char* pFilename)
{
  int nDim;
  int i, j;
  
  FILE* pFile = fopen(pFilename,"r");

  fscanf(pFile, "%d\n", &nDim);

  if(pMat->nD != nDim)
  {
    fclose(pFile);
    return 1;
  }

  for(i=0; i<nDim; i++)
  {
    for(j=0; j<nDim; j++)
    {
      fscanf(pFile, "%lf", &(pMat->m[i][j]));
    }
  }

  fclose(pFile);
 
  return 0;
}

//-------------------------------------------------------------------------------
// implement later
int MatrixTriangulation(SMatr* pMat, SMatr* pRes)
{
  return 0;
}

//-------------------------------------------------------------------------------
// implement later
int MatrixLU(SMatr* pMat, SMatr* pRes)
{
  return 0;
}