/* mat.c Simple matrix manipulation routines.
 * (C) Copyright 1999 by John Halleck
 * All rights reserved.
 */
/* Version of August 18th, 1999 */

/* Initial Test Version
 * 0.11 - John Halleck, September 6th
 *   added matrix subtract.
 * 0.10 - John Halleck, July 26th, 1999
 *   removed embed and extract, and moved them to a partitioning package.
 * 0.09 - John Halleck, July 24th, 1999
 *   Added embed and extract operations to support matrix partitioning.
 * 0.08 - John Halleck, July 17, 1999
 *   Stripped out most of the guts.  Pulled vector stuff.
 *   Only the stuff we use is left.  That should make
 *   it less work to make matching symmetric versions.
 * 0.07 - John Halleck, July 16, 1999
 *   Stripped some functions out.
 * 0.06 - John Halleck, June 30th, 1999
 *   Regularized form to match rest of package.
 * 0.05 - John Halleck, June 29th, 1999
 *   pulled various max functions.  Added copy transpose.
 * 0.04 - John Halleck, June 21st, 1999
 *   added rowmuls
 * 0.03 - John Halleck June 8th, 1999
 *   added some routines for interface completeness.
 * 0.02 - John Halleck June 5th, 1999
 *   added row vector vs col vector stuff.
 * 0.01 - John Halleck June 4, 1999
 *   Initial creation.
 */

#include "errors.h"
 /* Standard error codes. */

#include "mat.h"
 /* Matrix definitions. */

/* ========================= Set to a constant ============================== */

/* ------------- matzero ---------------- */
/* Make a zero matrix. */

error matzero (int rows, int cols, matrix result) {
   int size; double *out;
   if (!result)          return ERRnil;  /* We have to have a matrix */
   if (rows<1 || cols<1) return ERRsize; /* Of reasonable size       */

   size = rows*cols; out = (double *) result;
   while (size--) *out++ = 0.0;
   return NoError;
}

/* ------------- matidn ----------------- */

/* matrix identity */
error matidn (int size, matrix result)
{
  double *lmatrix;
  int index; long fullsize;

  if (!result)  return ERRnil;   /* Must have an array */
  if (size < 1) return ERRsize;  /* With an actual size */

  fullsize = size*size;

  index = fullsize;  lmatrix = (double *) result;
  while (index--)   *lmatrix++ = 0.0;

  index = size;      lmatrix = (double *) result;
  while (index--) { *lmatrix = 1.0; lmatrix += size + 1; }

  return NoError;
}

/* ========================= Copy =========================================== */

/* ------------- matcpy ----------------- */
/* Copy one matrix to another. */
error matcpy (int rows, int cols, matrix result, matrix source) {
    long size; double *final, *initial;

    if (!result || !source)   return ERRnil;   /* Must have arrays */
    if (rows < 1 || cols < 1) return ERRsize;  /* of reasonable size */

    size = rows*cols; final = (double *)result; initial = (double *) source;
    while (size--) { *final++ = *initial++; }

    return NoError;
}

/* ---------- Copy Transpose ------------ */
/* Copy a matrix to its transpose */
/* Source and destination must be distinct. */
/* rows and columns refer to destination. */
error matcpyt (int rows, int cols, matrix result, matrix source) {
    double *final, *initial;
    int i, j;

    if (!result || !source)   return ERRnil;   /* Must have arrays */
    if (rows < 1 || cols < 1) return ERRsize;  /* of reasonable size */
    if (source == result)     return ERRsame;  /* And they must be distinct */

    final = (double *)result; initial = (double *) source;
    for (i=0;i<rows;i++) for (j=0;j<cols;j++) {
      *final++ = *(initial + j*rows + i);
    }

    return NoError;
}

/* ======================== Multiply by scalar ============================== */

/* ------------------ matmuls ---------------- */
/* Multiply matrix by scalar. */
error matscl (int rows, int cols, matrix result, matrix given, double scalar) {
    long size; double *dst, *src;

    if (!result || !given)    return ERRnil;   /* Must have an array */
    if (rows < 1 || cols < 1) return ERRsize;  /* of reasonable size */

    dst = (double *) result;  src = (double *) given;
    size  = rows * cols;
    while (size--)  *dst++ =  scalar * *src++;

    if (scalar == 0.0) return ERRmeaning;
    return NoError;
}

/* ========================== Add =========================================== */
/* result = A + B */

/* ------------------ matadd ------------------ */
/* Add two matrices. */
error matadd (int rows, int cols, matrix result, matrix A, matrix B) {
    long size; double *final, *locala, *localb;

    if (!result || !A || !B)  return ERRnil;   /* Must have arrays */
    if (rows < 1 || cols < 1) return ERRsize;  /* of reasonable size */

    final = (double *) result; locala = (double *) A; localb = (double *) B;
    size  = rows * cols;
    while (size--) *final++ = *locala++ + *localb++; 
    return NoError;
}

/* result = A - B */
/* ----------------- matsub ------------------ */
error matsub (int rows, int cols, matrix result, matrix A, matrix B) {
    long size; double *final, *locala, *localb;

    if (!result || !A || !B)  return ERRnil;   /* Must have arrays */
    if (rows < 1 || cols < 1) return ERRsize;  /* of reasonable size */

    final = (double *) result; locala = (double *) A; localb = (double *) B;
    size  = rows * cols;
    while (size--) *final++ = *locala++ - *localb++;
    return NoError;
}

/*  ========================== Multiply ===================================== */

/* ------------------ matmul ----------------- */
/* Multiply  matrices.   result = matrix a multiplied by matrix b;  */
/* result MUST be a matrix distinct from either A or B. */
error matmul (int rrows, int rcols, matrix result,
            int isize, matrix A, matrix B)
/* result is rrows x rcols */
/* A      is rrows x isize */
/* B      is isize x rcols */

{
  double mterm; /* Term we are collecting */
  double *rowpoint, *colpoint; /* pointer to start of  current row and column */
  double *ipointa, *ipointb;   /* current locations for inner product. */
  double *final;

  int thisrow, thiscol, index;

  if (!result || !A || !B)        return ERRnil; /* Must have arrays */
  if (rrows < 1 || rcols < 1 || isize < 1)
                                  return ERRsize;  /* With real size */
  if (result == A || result == B) return ERRsame;  /* and result is distinct */

  final = (double *) result; rowpoint = (double *) A;
  for (thisrow = 0; thisrow < rrows; thisrow++) {
      colpoint = (double *) B;
      for (thiscol = 0; thiscol < rcols; thiscol++) {
          mterm = 0.0;
          ipointa = rowpoint; ipointb = colpoint;
          for (index=0; index < isize; index++) {
            mterm += *ipointa++ * *ipointb; ipointb += rcols;
          }
          *final++ = mterm;
          colpoint++;    /* On to next column */
      }
      rowpoint += isize; /* Pointer to first item of this row in A */
  }
  return NoError;
}

/* ------------------ matmult ----------------- */
/* Multiply  matrix by a transpose. */

error matmult (int rrows, int rcols, matrix result,
               int isize, matrix A, matrix B)
/* result is rrows x rcols    */
/* A      is rrows x isize    */
/* B      is rcols x isize <= */

{
  double mterm; /* Term we are collecting */
  double *rowpoint, *colpoint; /* pointer to start of  current row and column */
  double *ipointa,  *ipointb;  /* current locations for inner product. */
  double *final;

  int thisrow, thiscol, index;

  if (!result || !A || !B)        return ERRnil; /* Must have arrays */
  if (rrows < 1 || rcols < 1 || isize < 1)
                                  return ERRsize;  /* With real size */
  if (result == A || result == B) return ERRsame;  /* and result is distinct */

  final = (double *) result;
  rowpoint = (double *) A; colpoint = (double *) B;
  for (thisrow = 0; thisrow < rrows; thisrow++) {
      for (thiscol = 0; thiscol < rcols; thiscol++) {
          mterm = 0.0;
          ipointa = rowpoint; ipointb = colpoint;
          for (index=0; index < isize; index++) {
            mterm += *ipointa++ * *ipointb++;
          }
          *final++ = mterm;
          colpoint += isize;
      }
      rowpoint += isize; /* Pointer to first item of this row in A */
      colpoint = (double *) B;
  }
  return NoError;
}