/*

  Teem: Tools to process and visualize scientific data and images             .

  Copyright (C) 2013, 2012, 2011, 2010, 2009  University of Chicago

  Copyright (C) 2008, 2007, 2006, 2005  Gordon Kindlmann

  Copyright (C) 2004, 2003, 2002, 2001, 2000, 1999, 1998  University of Utah

  This library is free software; you can redistribute it and/or

  modify it under the terms of the GNU Lesser General Public License

  (LGPL) as published by the Free Software Foundation; either

  version 2.1 of the License, or (at your option) any later version.

  The terms of redistributing and/or modifying this software also

  include exceptions to the LGPL that facilitate static linking.

  This library is distributed in the hope that it will be useful,

  but WITHOUT ANY WARRANTY; without even the implied warranty of

  MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the GNU

  Lesser General Public License for more details.

  You should have received a copy of the GNU Lesser General Public License

  along with this library; if not, write to Free Software Foundation, Inc.,

  51 Franklin Street, Fifth Floor, Boston, MA  02110-1301  USA

*/

#include "ell.h"

int

ell_Nm_check(Nrrd *mat, int doNrrdCheck) {

  static const char me[]="ell_Nm_check";

    }

  } else {

    }

  }

  }

  }

/*

******** ell_Nm_tran

**

**     M             N

** N [trn]  <--  M [mat]

*/

int

ell_Nm_tran(Nrrd *ntrn, Nrrd *nmat) {

  static const char me[]="ell_Nm_tran";

  double *mat, *trn;

  size_t MM, NN, mm, nn;

  }

  }

  /*

  if (nrrdAxesSwap(ntrn, nmat, 0, 1)) {

    biffMovef(ELL, NRRD, "%s: trouble", me);

    return 1;

    }

  */

                        MM, NN)) {

  }

    }

  }

/*

******** ell_Nm_mul

**

** Currently, only useful for matrix-matrix multiplication

**

** matrix-matrix:      M       N

**                  L [A] . M [B]

*/

int

ell_Nm_mul(Nrrd *nAB, Nrrd *nA, Nrrd *nB) {

  static const char me[]="ell_Nm_mul";

  double *A, *B, *AB, tmp;

  size_t LL, MM, NN, ll, mm, nn;

  }

  }

  }

                        NN, LL)) {

  }

      tmp = 0;

      }

    }

  }

/*

** _ell_LU_decomp()

**

** in-place LU decomposition

*/

int

_ell_LU_decomp(double *aa, size_t *indx, size_t NN)  {

  static const char me[]="_ell_LU_decomp";

  int ret=0;

  size_t ii, imax=0, jj, kk;

  double big, sum, tmp;

  double *vv;

  }

  /* find vv[i]: max of abs of everything in column i */

    big = 0.0;

        big = tmp;

    }

      ret = 1; goto seeya;

    }

  }

    /* for aa[ii][jj] in lower triangle (below diagonal), subtract from

       aa[ii][jj] the dot product of all elements to its left with elements

       above it (starting at the top) */

      }

    }

    /* for aa[ii][jj] in upper triangle (including diagonal), subtract from

       aa[ii][jj] the dot product of all elements above it with elements to

       its left (starting from the left) */

    big = 0.0;

      }

      /* imax column is one in which abs(aa[i][j])/vv[i] */

        big = tmp;

        imax = ii;

    }

    /* unless we're on the imax column, swap this column the with imax column,

       and permute vv[] accordingly */

      /* could record parity # of permutes here */

      }

    }

    }

    /* divide everything right of a[jj][jj] by a[jj][jj] */

      }

    }

  }

 seeya:

/*

** _ell_LU_back_sub

**

** given the matrix and index array from _ellLUDecomp generated from

** some matrix M, solves for x in the linear equation Mx = b, and

** puts the result back into b

*/

void

_ell_LU_back_sub(double *aa, size_t *indx, double *bb, size_t NN) {

  size_t ii, jj;

  double sum;

  /* Forward substitution, with lower triangular matrix */

    }

  }

  /* Backward substitution, with upper triangular matrix */

    }

  }

  return;

/*

** _ell_inv

**

** Invert NNxNN matrix based on LU-decomposition

**

** The given matrix is copied, turned into its LU-decomposition, and

** then repeated backsubstitution is used to get successive columns of

** the inverse.

*/

int

_ell_inv(double *inv, double *_mat, size_t NN) {

  static const char me[]="_ell_inv";

  size_t ii, jj, *indx=NULL;

  double *col=NULL, *mat=NULL;

  int ret=0;

  }

  }

    /* set column jj of inv to result of backsub */

    }

  }

 seeya:

}

/*

******** ell_Nm_inv

**

** computes the inverse of given matrix in nmat, and puts the

** inverse in the (maybe allocated) ninv.  Does not touch the

** values in nmat.

*/

int

ell_Nm_inv(Nrrd *ninv, Nrrd *nmat) {

  static const char me[]="ell_Nm_inv";

  double *mat, *inv;

  size_t NN;

  }

    return 1;

  }

                        NN, NN)) {

  }

  }

/*

******** ell_Nm_pseudo_inv()

**

** determines the pseudoinverse of the given matrix M by using the formula

** P = (M^T * M)^(-1) * M^T

**

** I'll get an SVD-based solution working later, since that gives a more

** general solution

*/

int

ell_Nm_pseudo_inv(Nrrd *ninv, Nrrd *nA) {

  static const char me[]="ell_Nm_pseudo_inv";

  Nrrd *nAt, *nAtA, *nAtAi;

  int ret=0;

  }

  }

 seeya:

/*

******** ell_Nm_wght_pseudo_inv()

**

** determines a weighted least squares solution via

** P = (A^T * W * A)^(-1) * A^T * W

*/

int

ell_Nm_wght_pseudo_inv(Nrrd *ninv, Nrrd *nA, Nrrd *nW) {

  static const char me[]="ell_Nm_wght_pseudo_inv";

  Nrrd *nAt, *nAtW, *nAtWA, *nAtWAi;

  int ret=0;

  }

  }

 seeya:

}