esl_dmatrix.c

/* Linear algebra operations in double-precision matrices.
 * 
 * Implements ESL_DMATRIX (double-precision matrix) and 
 * ESL_PERMUTATION (permutation matrix) objects.
 * 
 * Table of contents:
 *   1. The ESL_DMATRIX object
 *   2. Debugging/validation code for ESL_DMATRIX
 *   3. Visualization tools
 *   4. The ESL_PERMUTATION object
 *   5. Debugging/validation code for ESL_PERMUTATION
 *   6. The rest of the dmatrix API
 *   7. Optional: Interoperability with GSL
 *   8. Optional: Interfaces to LAPACK
 *   9. Unit tests
 *  10. Test driver
 *  11. Examples
 *
 * To do:
 *   - eventually probably want additional matrix types
 *   - unit tests poor 
 */
#include "esl_config.h"

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

#include "easel.h"
#include "esl_vectorops.h"
#include "esl_dmatrix.h"


/*****************************************************************
 * 1. The ESL_DMATRIX object.
 *****************************************************************/

/* Function:  esl_dmatrix_Create()
 *
 * Purpose:   Creates a general <n> x <m> matrix (<n> rows, <m> 
 *            columns).
 *
 * Args:      <n> - number of rows;    $>= 1$
 *            <m> - number of columns; $>= 1$
 * 
 * Returns:   a pointer to a new <ESL_DMATRIX> object. Caller frees
 *            with <esl_dmatrix_Destroy()>.
 *
 * Throws:    <NULL> if an allocation failed.
 */
ESL_DMATRIX *
esl_dmatrix_Create(int n, int m)
{
  ESL_DMATRIX *A = NULL;
  int r;
  int status;

  ESL_ALLOC(A, sizeof(ESL_DMATRIX));
  A->mx = NULL;
  A->n  = n;
  A->m  = m;

  ESL_ALLOC(A->mx, sizeof(double *) * n);
  A->mx[0] = NULL;

  ESL_ALLOC(A->mx[0], sizeof(double) * n * m);
  for (r = 1; r < n; r++)
    A->mx[r] = A->mx[0] + r*m;

  A->type   = eslGENERAL;
  A->ncells = n * m; 
  return A;
  
 ERROR:
  esl_dmatrix_Destroy(A);
  return NULL;
}


/* Function:  esl_dmatrix_CreateUpper()
 * Incept:    SRE, Wed Feb 28 08:45:45 2007 [Janelia]
 *
 * Purpose:   Creates a packed upper triangular matrix of <n> rows and
 *            <n> columns. Caller may only access cells $i \leq j$.
 *            Cells $i > j$ are not stored and are implicitly 0.
 *            
 *            Not all matrix operations in Easel can work on packed
 *            upper triangular matrices.
 *
 * Returns:   a pointer to a new <ESL_DMATRIX> object of type
 *            <eslUPPER>. Caller frees with <esl_dmatrix_Destroy()>.
 *
 * Throws:    <NULL> if allocation fails.
 *
 * Xref:      J1/10
 */
ESL_DMATRIX *
esl_dmatrix_CreateUpper(int n)
{
  int status;
  ESL_DMATRIX *A = NULL;
  int r;			/* counter over rows */
  int nc;			/* cell counter */

  /* matrix structure allocation */
  ESL_ALLOC(A, sizeof(ESL_DMATRIX)); 
  A->mx = NULL;
  A->n  = n;
  A->m  = n;

  /* n row ptrs */
  ESL_ALLOC(A->mx, sizeof(double *) * n); 
  A->mx[0] = NULL;

  /* cell storage */
  ESL_ALLOC(A->mx[0], sizeof(double) * n * (n+1) / 2);
  
  /* row pointers set in a tricksy overlapping way, so
   * mx[i][j] access works normally but only i<=j are valid.
   * xref J1/10.
   */
  nc = n;  /* nc is the number of valid cells assigned to rows so far */
  for (r = 1; r < n; r++) {
    A->mx[r] = A->mx[0] + nc - r; /* -r overlaps this row w/ previous row */
    nc += n-r;
  }
  A->type   = eslUPPER;
  A->ncells = n * (n+1) / 2; 
  return A;

 ERROR:
  esl_dmatrix_Destroy(A);
  return NULL;
}


/* Function:  esl_dmatrix_Destroy()
 *            
 * Purpose:   Frees an <ESL_DMATRIX> object <A>.
 */
int
esl_dmatrix_Destroy(ESL_DMATRIX *A)
{
  if (A != NULL && A->mx != NULL && A->mx[0] != NULL) free(A->mx[0]);
  if (A != NULL && A->mx != NULL)                     free(A->mx);
  if (A != NULL)                                      free(A);
  return eslOK;
}


/* Function:  esl_dmatrix_Copy()
 *
 * Purpose:   Copies <src> matrix into <dest> matrix. <dest> must
 *            be allocated already by the caller.
 * 
 *            You may copy to a matrix of a different type, so long as
 *            the copy makes sense. If <dest> matrix is a packed type
 *            and <src> is not, the values that should be zeros must
 *            be zero in <src>, else the routine throws
 *            <eslEINCOMPAT>. If the <src> matrix is a packed type and
 *            <dest> is not, the values that are implicitly zeros are
 *            set to zeros in the <dest> matrix.
 *            
 * Returns:   <eslOK> on success.
 *
 * Throws:    <eslEINCOMPAT> if <src>, <dest> are different sizes,
 *            or if their types differ and <dest> cannot represent
 *            <src>.
 */
int
esl_dmatrix_Copy(const ESL_DMATRIX *src, ESL_DMATRIX *dest)
{
  int i,j;

  if (dest->n != src->n || dest->m != src->m)
    ESL_EXCEPTION(eslEINCOMPAT, "matrices of different size");

  if (src->type == dest->type)   /* simple case. */
    memcpy(dest->mx[0], src->mx[0], src->ncells * sizeof(double));

  else if (src->type == eslGENERAL && dest->type == eslUPPER)		
    {
      for (i = 1; i < src->n; i++)
	for (j = 0; j < i; j++)
	  if (src->mx[i][j] != 0.) 
	    ESL_EXCEPTION(eslEINCOMPAT, "general matrix isn't upper triangular, can't be copied/packed");
      for (i = 0; i < src->n; i++)
	for (j = i; j < src->m; j++)
	  dest->mx[i][j] = src->mx[i][j];
    }
  
  else if (src->type == eslUPPER && dest->type == eslGENERAL)		
    {
      for (i = 1; i < src->n; i++)
	for (j = 0; j < i; j++)
	  dest->mx[i][j] = 0.;
      for (i = 0; i < src->n; i++)
	for (j = i; j < src->m; j++)
	  dest->mx[i][j] = src->mx[i][j];      
    }

  return eslOK;
}


/* Function:  esl_dmatrix_Clone()
 * Incept:    SRE, Tue May  2 14:38:45 2006 [St. Louis]
 *
 * Purpose:   Duplicates matrix <A>, making a copy in newly
 *            allocated space.
 *
 * Returns:   a pointer to the copy. Caller frees with 
 *            <esl_dmatrix_Destroy()>.
 *
 * Throws:    <NULL> on allocation failure.
 */
ESL_DMATRIX *
esl_dmatrix_Clone(const ESL_DMATRIX *A)
{
  ESL_DMATRIX *new;

  switch (A->type) {
  case eslUPPER:             if ( (new = esl_dmatrix_CreateUpper(A->n))  == NULL) return NULL; break;
  default: case eslGENERAL:  if ( (new = esl_dmatrix_Create(A->n, A->m)) == NULL) return NULL; break;
  }
  esl_dmatrix_Copy(A, new);
  return new;
}


/* Function:  esl_dmatrix_Compare()
 *
 * Purpose:   Compares matrix <A> to matrix <B> element by element,
 *            using <esl_DCompare()> on each cognate element pair,
 *            with relative equality defined by a fractional tolerance
 *            <tol>.  If all elements are equal, return <eslOK>; if
 *            any elements differ, return <eslFAIL>.
 *            
 *            <A> and <B> may be of different types; for example,
 *            a packed upper triangular matrix A is compared to
 *            a general matrix B by assuming <A->mx[i][j] = 0.> for
 *            all $i>j$.
 */
int
esl_dmatrix_Compare(const ESL_DMATRIX *A, const ESL_DMATRIX *B, double tol)
{
  int i,j,c;
  double x1,x2;

  if (A->n != B->n) return eslFAIL;
  if (A->m != B->m) return eslFAIL;

  if (A->type == B->type) 
    {  /* simple case. */
      for (c = 0; c < A->ncells; c++) /* can deal w/ packed or unpacked storage */
	if (esl_DCompare(A->mx[0][c], B->mx[0][c], tol) == eslFAIL) return eslFAIL;
    }
  else 
    { /* comparing matrices of different types */
      for (i = 0; i < A->n; i++)
	for (j = 0; j < A->m; j++)
	  {
	    if (A->type == eslUPPER && i > j) x1 = 0.;
	    else                              x1 = A->mx[i][j];

	    if (B->type == eslUPPER && i > j) x2 = 0.;
	    else                              x2 = B->mx[i][j];

	    if (esl_DCompare(x1, x2, tol) == eslFAIL) return eslFAIL;
	  }
    }
  return eslOK;
}


/* Function:  esl_dmatrix_CompareAbs()
 *
 * Purpose:   Compares matrix <A> to matrix <B> element by element,
 *            using <esl_DCompareAbs()> on each cognate element pair,
 *            with absolute equality defined by a absolute difference tolerance
 *            <tol>.  If all elements are equal, return <eslOK>; if
 *            any elements differ, return <eslFAIL>.
 *            
 *            <A> and <B> may be of different types; for example,
 *            a packed upper triangular matrix A is compared to
 *            a general matrix B by assuming <A->mx[i][j] = 0.> for
 *            all $i>j$.
 */
int
esl_dmatrix_CompareAbs(const ESL_DMATRIX *A, const ESL_DMATRIX *B, double tol)
{
  int i,j,c;
  double x1,x2;

  if (A->n != B->n) return eslFAIL;
  if (A->m != B->m) return eslFAIL;

  if (A->type == B->type) 
    {  /* simple case. */
      for (c = 0; c < A->ncells; c++) /* can deal w/ packed or unpacked storage */
	if (esl_DCompareAbs(A->mx[0][c], B->mx[0][c], tol) == eslFAIL) return eslFAIL;
    }
  else 
    { /* comparing matrices of different types */
      for (i = 0; i < A->n; i++)
	for (j = 0; j < A->m; j++)
	  {
	    if (A->type == eslUPPER && i > j) x1 = 0.;
	    else                              x1 = A->mx[i][j];

	    if (B->type == eslUPPER && i > j) x2 = 0.;
	    else                              x2 = B->mx[i][j];

	    if (esl_DCompareAbs(x1, x2, tol) == eslFAIL) return eslFAIL;
	  }
    }
  return eslOK;
}


/* Function:  esl_dmatrix_Set()
 *
 * Purpose:   Set all elements $a_{ij}$ in matrix <A> to <x>,
 *            and returns <eslOK>.
 */
int
esl_dmatrix_Set(ESL_DMATRIX *A, double x)
{
  int i;
  for (i = 0; i < A->ncells; i++) A->mx[0][i] = x;
  return eslOK;
}


/* Function:  esl_dmatrix_SetZero()
 *
 * Purpose:   Sets all elements $a_{ij}$ in matrix <A> to 0,
 *            and returns <eslOK>.
 */
int
esl_dmatrix_SetZero(ESL_DMATRIX *A)
{
  int i;
  for (i = 0; i < A->ncells; i++) A->mx[0][i] = 0.;
  return eslOK;
}
  

/* Function:  esl_dmatrix_SetIdentity()
 *
 * Purpose:   Given a square matrix <A>, sets all diagonal elements 
 *            $a_{ii}$ to 1, and all off-diagonal elements $a_{ij},
 *            j \ne i$ to 0. Returns <eslOK> on success.
 *
 * Throws:    <eslEINVAL> if the matrix isn't square.
 */
int
esl_dmatrix_SetIdentity(ESL_DMATRIX *A)
{
  int i;
  
  if (A->n != A->m) ESL_EXCEPTION(eslEINVAL, "matrix isn't square");
  esl_dmatrix_SetZero(A);
  for (i = 0; i < A->n; i++) A->mx[i][i] = 1.;
  return eslOK;
}
  
/*****************************************************************
 * 2. Debugging, validation code
 *****************************************************************/

/* Function:  esl_dmatrix_Dump()
 * Incept:    SRE, Mon Nov 29 19:21:20 2004 [St. Louis]
 *
 * Purpose:   Given a matrix <A>, dump it to output stream <ofp> in human-readable
 *            format.
 * 
 *            If <rowlabel> or <collabel> are non-NULL, they specify a
 *            string of single-character labels to put on the rows and
 *            columns, respectively. (For example, these might be a
 *            sequence alphabet for a 4x4 or 20x20 rate matrix or
 *            substitution matrix.)  Numbers <1..ncols> or <1..nrows> are
 *            used if <collabel> or <rowlabel> are passed as <NULL>.
 *
 * Args:      ofp      -  output file pointer; stdout, for example.
 *            A        -  matrix to dump.
 *            rowlabel -  optional: NULL, or character labels for rows
 *            collabel -  optional: NULL, or character labels for cols
 *
 * Returns:   <eslOK> on success.
 */
int
esl_dmatrix_Dump(FILE *ofp, const ESL_DMATRIX *A, const char *rowlabel, const char *collabel)
{
  int a,b;

  fprintf(ofp, "     ");
  if (collabel != NULL) 
    for (b = 0; b < A->m; b++) fprintf(ofp, "       %c ", collabel[b]);
  else
    for (b = 0; b < A->m; b++) fprintf(ofp, "%8d ", b+1);
  fprintf(ofp, "\n");

  for (a = 0; a < A->n; a++) {
    if (rowlabel != NULL)      fprintf(ofp, "    %c ", rowlabel[a]);
    else                       fprintf(ofp, "%5d ",    a+1);

    for (b = 0; b < A->m; b++) {
      switch (A->type) {
      case eslUPPER:
	if (a > b) 	fprintf(ofp, "%8s ", "");
	else            fprintf(ofp, "%8.4f ", A->mx[a][b]); 
	break;

       default: case eslGENERAL:
	fprintf(ofp, "%8.4f ", A->mx[a][b]); 
	break;
      }
    }
    fprintf(ofp, "\n");
  }
  return eslOK;
}


/*****************************************************************
 * 3. Visualization tools
 *****************************************************************/

/* Function:  esl_dmatrix_PlotHeatMap()
 * Synopsis:  Export a heat map visualization, in PostScript
 *
 * Purpose:   Export a heat map visualization of the matrix in <D>
 *            to open stream <fp>, in PostScript format. 
 *            
 *            All values between <min> and <max> in <D> are rescaled
 *            linearly and assigned to shades. Values below <min>
 *            are assigned to the lowest shade; values above <max>, to
 *            the highest shade.
 *
 *            The plot is hardcoded to be a full US 8x11.5" page,
 *            with at least a 20pt margin.
 *
 *            Several color schemes are enumerated in the code
 *            but all but one is commented out. The currently enabled
 *            scheme is a 10-class scheme consisting of the 9-class
 *            Reds from colorbrewer2.org plus a blue background class.
 *          
 * Note:      Binning rules basically follow same convention as
 *            esl_histogram. nb = xmax-xmin/w, so w = xmax-xmin/nb; 
 *            picking bin is (int) ceil((x - xmin)/w) - 1. (xref
 *            esl_histogram_Score2Bin()). This makes bin b contain
 *            values bw+min < x <= (b+1)w+min. (Which means that 
 *            min itself falls in bin -1, whoops - but we catch
 *            all bin<0 and bin>=nshades and put them in the extremes.
 *
 * Returns:   <eslOK> on success.
 *
 * Throws:    (no abnormal error conditions)
 */
int
esl_dmatrix_PlotHeatMap(FILE *fp, ESL_DMATRIX *D, double min, double max)
{
#if 0
  /*
   * This color scheme roughly follows Tufte, Envisioning Information,
   * p.91, where he shows a beautiful bathymetric chart. The CMYK
   * values conjoin two recommendations from ColorBrewer (Cindy Brewer
   * and Mark Harrower, colorbrewer2.org), specifically the 9-class
   * sequential2 Blues and 9-class sequential YlOrBr.
   */
  int    nshades   = 18;
  double cyan[]    = { 1.00, 1.00, 0.90, 0.75, 0.57, 0.38, 0.24, 0.13, 0.03,
                       0.00, 0.00, 0.00, 0.00, 0.00, 0.07, 0.20, 0.40, 0.60};
  double magenta[] = { 0.55, 0.45, 0.34, 0.22, 0.14, 0.08, 0.06, 0.03, 0.01,
                       0.00, 0.03, 0.11, 0.23, 0.40, 0.55, 0.67, 0.75, 0.80};
  double yellow[]  = { 0.00, 0.00, 0.00, 0.00, 0.00, 0.00, 0.00, 0.00, 0.00,
                       0.10, 0.25, 0.40, 0.65, 0.80, 0.90, 1.00, 1.00, 1.00};
  double black[]   = { 0.30, 0.07, 0.00, 0.00, 0.00, 0.00, 0.00, 0.00, 0.00,
                       0.00, 0.00, 0.00, 0.00, 0.00, 0.00, 0.00, 0.00, 0.00};
#endif
#if 0
  /* colorbrewer2.org 5-class YlOrBr scheme: sequential, multihue, 5-class, CMYK */
  int    nshades   = 5;
  double cyan[]    = { 0.00, 0.00, 0.00, 0.15, 0.40 };
  double magenta[] = { 0.00, 0.15, 0.40, 0.60, 0.75 };
  double yellow[]  = { 0.17, 0.40, 0.80, 0.95, 1.00 };
  double black[]   = { 0,    0,    0,    0,    0    };
#endif
#if 0
  /* colorbrewer2.org 9-class YlOrBr scheme, +zero class */
  int    nshades   = 10;
  double cyan[]    = { 0.00, 0.00, 0.00, 0.00, 0.00, 0.00, 0.07, 0.20, 0.40, 0.60 };
  double magenta[] = { 0.00, 0.00, 0.03, 0.11, 0.23, 0.40, 0.55, 0.67, 0.75, 0.80 };
  double yellow[]  = { 0.00, 0.10, 0.25, 0.40, 0.65, 0.80, 0.90, 1.00, 1.00, 1.00 };
  double black[]   = { 0.05, 0.00, 0.00, 0.00, 0.00, 0.00, 0.00, 0.00, 0.00, 0.00 };
#endif
  /* colorbrewer2.org 9-class Reds + zero class as dim blue */
  int    nshades   = 10;
  double cyan[]    = { 0.30, 0.00, 0.00, 0.00, 0.00, 0.00, 0.05, 0.20, 0.35, 0.60 };
  double magenta[] = { 0.03, 0.04, 0.12, 0.27, 0.43, 0.59, 0.77, 0.90, 0.95, 1.00 };
  double yellow[]  = { 0.00, 0.04, 0.12, 0.27, 0.43, 0.59, 0.72, 0.80, 0.85, 0.90 };
  double black[]   = { 0.00, 0.00, 0.00, 0.00, 0.00, 0.00, 0.00, 0.00, 0.00, 0.00 };

  int    pageheight = 792;
  int    pagewidth  = 612;
  double w;                     
  int    i,j;
  int    bin;
  float  boxsize;               /* box size in points */
  float  xcoord, ycoord;        /* postscript coords in points */
  float  leftmargin;
  float  bottommargin;

  /* Set some defaults that might become arguments later.
   */
  leftmargin  = 20.;
  bottommargin = 20.;

  /* Determine some working parameters 
   */
  w = (max-min) / (double) nshades; /* w = bin size for assigning values->colors*/
  boxsize = ESL_MIN( (pageheight - (bottommargin * 2.)) / (float) D->n, 
                     (pagewidth - (leftmargin * 2.))   / (float) D->m);
  
  /* or start from j=i, to do diagonals */
  for (i = 0; i < D->n; i++)
    for (j = 0; j < D->m; j++)
      {
        xcoord = (float) j * boxsize + leftmargin;
        ycoord = (float) (D->n-i+1) * boxsize + bottommargin;

        if      (D->mx[i][j] == -eslINFINITY) bin = 0;
        else if (D->mx[i][j] ==  eslINFINITY) bin = nshades-1;
        else {
          bin    = (int) ceil((D->mx[i][j] - min) / w) - 1;
          if (bin < 0)        bin = 0;
          if (bin >= nshades) bin = nshades-1;
        }

        fprintf(fp, "newpath\n");
        fprintf(fp, "  %.2f %.2f moveto\n", xcoord, ycoord);
        fprintf(fp, "  0  %.2f rlineto\n", boxsize);
        fprintf(fp, "  %.2f 0  rlineto\n", boxsize);
        fprintf(fp, "  0 -%.2f rlineto\n", boxsize);
        fprintf(fp, "  closepath\n");
        fprintf(fp, " %.2f %.2f %.2f %.2f setcmykcolor\n",
                cyan[bin], magenta[bin], yellow[bin], black[bin]);
        fprintf(fp, "  fill\n");
      }
  fprintf(fp, "showpage\n");
  return eslOK;
}



/*****************************************************************
 * 4. The ESL_PERMUTATION object.
 *****************************************************************/

/* Function:  esl_permutation_Create()
 *
 * Purpose:   Creates a new permutation "matrix" of size <n> for
 *            permuting <n> x <n> square matrices; returns a 
 *            pointer to it.
 *
 *            A permutation matrix consists of 1's and 0's such that
 *            any given row or column contains only one 1. We store it
 *            more efficiently as a vector; each value $p_i$
 *            represents the column $j$ that has the 1. Thus, on
 *            initialization, $p_i = i$ for all $i = 0..n-1$.
 *
 * Returns:   a pointer to a new <ESL_PERMUTATION> object. Free with 
 *            <esl_permutation_Destroy()>.
 *
 * Throws:    <NULL> if allocation fails.
 */
ESL_PERMUTATION *
esl_permutation_Create(int n)
{
  int status;
  ESL_PERMUTATION *P = NULL;

  ESL_DASSERT1(( n > 0 ));

  ESL_ALLOC(P, sizeof(ESL_PERMUTATION));
  P->pi = NULL;
  P->n  = n;
  ESL_ALLOC(P->pi, sizeof(int) * n);

  esl_permutation_Reuse(P);	/* initialize it */
  return P;

 ERROR:
  esl_permutation_Destroy(P);
  return NULL;
}
  
/* Function:  esl_permutation_Destroy()
 *
 * Purpose:   Frees an <ESL_PERMUTATION> object <P>.
 */
int
esl_permutation_Destroy(ESL_PERMUTATION *P)
{
  if (P != NULL && P->pi != NULL) free(P->pi);
  if (P != NULL)                  free(P);
  return eslOK;
}

/* Function:  esl_permutation_Reuse()
 *
 * Purpose:   Resets a permutation matrix <P> to
 *            $p_i = i$ for all $i = 0..n-1$.
 *            
 * Returns:   <eslOK> on success.           
 */
int
esl_permutation_Reuse(ESL_PERMUTATION *P)
{
  int i;
  for (i = 0; i < P->n; i++)
    P->pi[i] = i;
  return eslOK;
}


/*****************************************************************
 * 5. Debugging/validation for ESL_PERMUTATION.
 *****************************************************************/

/* Function:  esl_permutation_Dump()
 *
 * Purpose:   Given a permutation matrix <P>, dump it to output stream <ofp>
 *            in human-readable format.
 *            
 *            If <rowlabel> or <collabel> are non-NULL, they represent
 *            single-character labels to put on the rows and columns,
 *            respectively. (For example, these might be a sequence
 *            alphabet for a 4x4 or 20x20 rate matrix or substitution
 *            matrix.)  Numbers 1..ncols or 1..nrows are used if
 *            <collabel> or <rowlabel> are NULL.
 *
 * Args:      ofp      - output file pointer; stdout, for example
 *            P        - permutation matrix to dump
 *            rowlabel - optional: NULL, or character labels for rows
 *            collabel - optional: NULL, or character labels for cols
 *
 * Returns:   <eslOK> on success.
 */
int 
esl_permutation_Dump(FILE *ofp, const ESL_PERMUTATION *P, const char *rowlabel, const char *collabel)
{
  int i,j;

  fprintf(ofp, "    ");
  if (collabel != NULL)
    for (j = 0; j < P->n; j++) fprintf(ofp, "  %c ", collabel[j]);
  else
    for (j = 0; j < P->n; j++) fprintf(ofp, "%3d ", j+1);
  fprintf(ofp, "\n");

  for (i = 0; i < P->n; i++) {
    if (rowlabel != NULL) fprintf(ofp, "  %c ", rowlabel[i]);
    else                  fprintf(ofp, "%3d ", i+1);

    for (j = 0; j < P->n; j++)
      fprintf(ofp, "%3d ", (j == P->pi[i]) ? 1 : 0);
    fprintf(ofp, "\n");
  }
  return eslOK;
}

/*****************************************************************
 * 6. The rest of the dmatrix API.
 *****************************************************************/



/* Function:  esl_dmx_Max()
 * Incept:    SRE, Thu Mar  1 14:46:48 2007 [Janelia]
 *
 * Purpose:   Returns the maximum value of all the elements $a_{ij}$ in matrix <A>.
 */
double
esl_dmx_Max(const ESL_DMATRIX *A)
{
  int    i;
  double best;

  best = A->mx[0][0];
  for (i = 0; i < A->ncells; i++)
    if (A->mx[0][i] > best) best = A->mx[0][i];
  return best;
}

/* Function:  esl_dmx_Min()
 * Incept:    SRE, Thu Mar  1 14:49:29 2007 [Janelia]
 *
 * Purpose:   Returns the minimum value of all the elements $a_{ij}$ in matrix <A>.
 */
double
esl_dmx_Min(const ESL_DMATRIX *A)
{
  int    i;
  double best;

  best = A->mx[0][0];
  for (i = 0; i < A->ncells; i++)
    if (A->mx[0][i] < best) best = A->mx[0][i];
  return best;
}


/* Function:  esl_dmx_MinMax()
 * Incept:    SRE, Wed Mar 14 16:58:03 2007 [Janelia]
 *
 * Purpose:   Finds the maximum and minimum values of the
 *            elements $a_{ij}$ in matrix <A>, and returns
 *            them in <ret_min> and <ret_max>.
 *            
 * Returns:   <eslOK> on success.            
 *            
 */
int
esl_dmx_MinMax(const ESL_DMATRIX *A, double *ret_min, double *ret_max)
{
  double min, max;
  int i;

  min = max = A->mx[0][0];
  for (i = 0; i < A->ncells; i++) {
    if (A->mx[0][i] < min) min = A->mx[0][i];
    if (A->mx[0][i] > max) max = A->mx[0][i];
  }
  *ret_min = min;
  *ret_max = max;
  return eslOK;
}



/* Function:  esl_dmx_Sum()
 * Incept:    SRE, Thu Mar  1 16:45:16 2007
 *
 * Purpose:   Returns the scalar sum of all the elements $a_{ij}$ in matrix <A>,
 *            $\sum_{ij} a_{ij}$.
 */
double
esl_dmx_Sum(const ESL_DMATRIX *A)
{
  int    i;
  double sum = 0.;

  for (i = 0; i < A->ncells; i++)
    sum += A->mx[0][i];
  return sum;
}


/* Function:  esl_dmx_FrobeniusNorm()
 * Incept:    SRE, Thu Mar 15 17:59:35 2007 [Janelia]
 *
 * Purpose:   Calculates the Frobenius norm of a matrix, which
 *            is the element-wise equivalant of a 
 *            Euclidean vector norm: 
 *               $ = \sqrt(\sum a_{ij}^2)$
 *
 * Args:      A         - matrix
 *            ret_fnorm - Frobenius norm.
 * 
 * Returns:   <eslOK> on success, and the Frobenius norm
 *            is in <ret_fnorm>.
 */
int
esl_dmx_FrobeniusNorm(const ESL_DMATRIX *A, double *ret_fnorm)
{
  double F = 0.;
  int    i;

  for (i = 0; i < A->ncells; i++)
    F += A->mx[0][i] * A->mx[0][i];
  *ret_fnorm = sqrt(F);
  return eslOK;
}


/* Function: esl_dmx_Multiply()
 * 
 * Purpose:  Matrix multiplication: calculate <AB>, store result in <C>.
 *           <A> is $n times m$; <B> is $m \times p$; <C> is $n \times p$.
 *           Matrix <C> must be allocated appropriately by the caller.
 *
 *           Not supported for anything but general (<eslGENERAL>)
 *           matrix type, at present.
 *           
 * Throws:   <eslEINVAL> if matrices don't have compatible dimensions,
 *           or if any of them isn't a general (<eslGENERAL>) matrix.
 */
int
esl_dmx_Multiply(const ESL_DMATRIX *A, const ESL_DMATRIX *B, ESL_DMATRIX *C)
{
  int i, j, k;

  if (A->m    != B->n)       ESL_EXCEPTION(eslEINVAL, "can't multiply A,B");
  if (A->n    != C->n)       ESL_EXCEPTION(eslEINVAL, "A,C # of rows not equal");
  if (B->m    != C->m)       ESL_EXCEPTION(eslEINVAL, "B,C # of cols not equal");
  if (A->type != eslGENERAL) ESL_EXCEPTION(eslEINVAL, "A isn't of type eslGENERAL");
  if (B->type != eslGENERAL) ESL_EXCEPTION(eslEINVAL, "B isn't of type eslGENERAL");
  if (C->type != eslGENERAL) ESL_EXCEPTION(eslEINVAL, "B isn't of type eslGENERAL");

  /* i,k,j order should optimize stride, relative to a more textbook
   * order for the indices
   */
  esl_dmatrix_SetZero(C);
  for (i = 0; i < A->n; i++)
    for (k = 0; k < A->m; k++)
      for (j = 0; j < B->m; j++)
	C->mx[i][j] += A->mx[i][k] * B->mx[k][j];

  return eslOK;
}


/*::cexcerpt::function_comment_example::begin::*/
/* Function:  esl_dmx_Exp()
 * Synopsis:  Calculates matrix exponential $\mathbf{P} = e^{t\mathbf{Q}}$.
 * Incept:    SRE, Thu Mar  8 18:41:38 2007 [Janelia]
 *
 * Purpose:   Calculates the matrix exponential $\mathbf{P} = e^{t\mathbf{Q}}$,
 *            using a scaling and squaring algorithm with
 *            the Taylor series approximation \citep{MolerVanLoan03}.
 *                              
 *            <Q> must be a square matrix of type <eslGENERAL>.
 *            Caller provides an allocated <P> matrix of the same size and type as <Q>.
 *            
 *            A typical use of this function is to calculate a
 *            conditional substitution probability matrix $\mathbf{P}$
 *            (whose elements $P_{xy}$ are conditional substitution
 *            probabilities $\mathrm{Prob}(y \mid x, t)$ from time $t$
 *            and instantaneous rate matrix $\mathbf{Q}$.
 *
 * Args:      Q  - matrix to exponentiate (an instantaneous rate matrix)
 *            t  - time units
 *            P  - RESULT: $e^{tQ}$.
 *
 * Returns:   <eslOK> on success.
 *
 * Throws:    <eslEMEM> on allocation error.
 *
 * Xref:      J1/19.
 */
int
esl_dmx_Exp(const ESL_DMATRIX *Q, double t, ESL_DMATRIX *P)
{
/*::cexcerpt::function_comment_example::end::*/
  ESL_DMATRIX *Qz   = NULL;	/* Q/2^z rescaled matrix*/
  ESL_DMATRIX *Qpow = NULL;	/* keeps running product Q^k */
  ESL_DMATRIX *C    = NULL;	/* tmp storage for matrix multiply result */
  double factor     = 1.0;
  double fnorm;
  int    z;
  double zfac;
  int    k;
  int    status;
    
  /* Contract checks  */
  if (Q->type != eslGENERAL) ESL_EXCEPTION(eslEINVAL, "Q isn't general");
  if (Q->n    != Q->m)       ESL_EXCEPTION(eslEINVAL, "Q isn't square");
  if (P->type != Q->type)    ESL_EXCEPTION(eslEINVAL, "P isn't of same type as Q");
  if (P->n    != P->m)       ESL_EXCEPTION(eslEINVAL, "P isn't square");
  if (P->n    != Q->n)       ESL_EXCEPTION(eslEINVAL, "P isn't same size as Q");

  /* Allocation of working space */
  if ((Qz   = esl_dmatrix_Create(Q->n, Q->n)) == NULL) { status = eslEMEM; goto ERROR; }
  if ((Qpow = esl_dmatrix_Create(Q->n, Q->n)) == NULL) { status = eslEMEM; goto ERROR; }
  if ((C    = esl_dmatrix_Create(Q->n, Q->n)) == NULL) { status = eslEMEM; goto ERROR; }
  
  /* Figure out how much to scale the matrix down by.  This is not
   * magical; we're just knocking its magnitude down in an ad hoc way.
   */
  esl_dmx_FrobeniusNorm(Q, &fnorm);
  zfac = 1.;
  z    = 0;
  while (t*fnorm*zfac > 0.1) { zfac /= 2.; z++; }

  /* Make a scaled-down copy of Q in Qz. 
   */ 
  esl_dmatrix_Copy(Q, Qz);       
  esl_dmx_Scale(Qz, zfac);

  /* Calculate e^{t Q_z} by the Taylor, to complete convergence. */
  esl_dmatrix_SetIdentity(P);
  esl_dmatrix_Copy(Qz, Qpow);       /* Qpow is now Qz^1 */
  for (k = 1; k < 100; k++)
    {
      factor *= t/k;
      esl_dmatrix_Copy(P, C);	            /* C now holds the previous P */
      esl_dmx_AddScale(P, factor, Qpow);    /* P += factor*Qpow */
      if (esl_dmatrix_Compare(C, P, 0.) == eslOK) break;

      esl_dmx_Multiply(Qpow, Qz, C);        /* C = Q^{k+1} */
      esl_dmatrix_Copy(C, Qpow);            /* Qpow = C = Q^{k+1} */
    }

  /* Now square it back up: e^{tQ} = [e^{tQ_z}]^{2^z} */
  while (z--) {
    esl_dmx_Multiply(P, P, C);
    esl_dmatrix_Copy(C, P);
  }

  esl_dmatrix_Destroy(Qz);
  esl_dmatrix_Destroy(Qpow);
  esl_dmatrix_Destroy(C);
  return eslOK;

 ERROR:
  if (Qz   != NULL) esl_dmatrix_Destroy(Qz);
  if (Qpow != NULL) esl_dmatrix_Destroy(Qpow);
  if (C    != NULL) esl_dmatrix_Destroy(C);
  return status;
}


/* Function:  esl_dmx_Transpose()
 *
 * Purpose:   Transpose a square matrix <A> in place.
 *
 *            <A> must be a general (<eslGENERAL>) matrix type.
 *
 * Throws:    <eslEINVAL> if <A> isn't square, or if it isn't
 *            of type <eslGENERAL>.
 */
int
esl_dmx_Transpose(ESL_DMATRIX *A)
{
  int    i,j;
  double swap;

  if (A->n    != A->m)       ESL_EXCEPTION(eslEINVAL, "matrix isn't square");
  if (A->type != eslGENERAL) ESL_EXCEPTION(eslEINVAL, "A isn't of type eslGENERAL");

  for (i = 0; i < A->n; i++)
    for (j = i+1; j < A->m; j++)
      { swap = A->mx[i][j]; A->mx[i][j] = A->mx[j][i]; A->mx[j][i] = swap; }
  return eslOK;
}


/* Function:  esl_dmx_Add()
 *
 * Purpose:   <A = A+B>; adds matrix <B> to matrix <A> and leaves result
 *            in matrix <A>.
 *
 *            <A> and <B> may be of any type. However, if <A> is a
 *            packed upper triangular matrix (type
 *            <eslUPPER>), all values $i>j$ in <B> must be
 *            zero (i.e. <B> must also be upper triangular, though
 *            not necessarily packed upper triangular).
 *
 * Throws:    <eslEINVAL> if matrices aren't the same dimensions, or
 *            if <A> is <eslUPPER> and any cell $i>j$ in
 *            <B> is nonzero.
 */
int
esl_dmx_Add(ESL_DMATRIX *A, const ESL_DMATRIX *B)
{
  int    i,j;
  
  if (A->n    != B->n)              ESL_EXCEPTION(eslEINVAL, "matrices of different size");
  if (A->m    != B->m)              ESL_EXCEPTION(eslEINVAL, "matrices of different size");

  if (A->type == B->type)	/* in this case, can just add cell by cell */
    {
      for (i = 0; i < A->ncells; i++)
	A->mx[0][i] += B->mx[0][i];
    }
  else if (A->type == eslUPPER || B->type == eslUPPER)
    {
      /* Logic is: if either matrix is upper triangular, then the operation is
       * to add upper triangles only. If we try to add a general matrix <B>
       * to packed UT <A>, make sure all lower triangle entries in <B> are zero.
       */
      if (B->type != eslUPPER) {
	for (i = 1; i < A->n; i++)
	  for (j = 0; j < i; j++)
	    if (B->mx[i][j] != 0.) ESL_EXCEPTION(eslEINVAL, "<B> has nonzero cells in lower triangle");
      }
      for (i = 0; i < A->n; i++)
	for (j = i; j < A->m; j++)
	  A->mx[i][j] += B->mx[i][j];
    }
  return eslOK;
}

/* Function:  esl_dmx_Scale()
 *
 * Purpose:   Calculates <A = kA>: multiply matrix <A> by scalar
 *            <k> and leave answer in <A>.
 */
int 
esl_dmx_Scale(ESL_DMATRIX *A, double k)
{
  int i;

  for (i = 0; i < A->ncells; i++)  A->mx[0][i] *=  k;
  return eslOK;
}


/* Function:  esl_dmx_AddScale()
 * 
 * Purpose:   Calculates <A + kB>, leaves answer in <A>.
 * 
 *            Only defined for matrices of the same type (<eslGENERAL>
 *            or <eslUPPER>).
 * 
 * Throws:    <eslEINVAL> if matrices aren't the same dimensions, or
 *            of different types.
 */
int
esl_dmx_AddScale(ESL_DMATRIX *A, double k, const ESL_DMATRIX *B)
{
  int i;

  if (A->n    != B->n)    ESL_EXCEPTION(eslEINVAL, "matrices of different size");
  if (A->m    != B->m)    ESL_EXCEPTION(eslEINVAL, "matrices of different size");
  if (A->type != B->type) ESL_EXCEPTION(eslEINVAL, "matrices of different type");

  for (i = 0; i < A->ncells; i++) A->mx[0][i] +=  k * B->mx[0][i];
  return eslOK;
}


/* Function:  esl_dmx_Permute_PA()
 *
 * Purpose:   Computes <B = PA>: do a row-wise permutation of a square
 *            matrix <A>, using the permutation matrix <P>, and put
 *            the result in a square matrix <B> that the caller has
 *            allocated.
 *
 * Throws:    <eslEINVAL> if <A>, <B>, <P> do not have compatible dimensions,
 *            or if <A> or <B> is not of type <eslGENERAL>.
 */
int
esl_dmx_Permute_PA(const ESL_PERMUTATION *P, const ESL_DMATRIX *A, ESL_DMATRIX *B)
{
  int i,ip,j;

  if (A->n    != P->n)       ESL_EXCEPTION(eslEINVAL, "matrix dimensions not compatible");
  if (A->n    != B->n)       ESL_EXCEPTION(eslEINVAL, "matrix dimensions not compatible");
  if (A->n    != A->m)       ESL_EXCEPTION(eslEINVAL, "matrix dimensions not compatible");
  if (B->n    != B->m)       ESL_EXCEPTION(eslEINVAL, "matrix dimensions not compatible");
  if (A->type != eslGENERAL) ESL_EXCEPTION(eslEINVAL, "matrix A not of type eslGENERAL");
  if (B->type != eslGENERAL) ESL_EXCEPTION(eslEINVAL, "matrix B not of type eslGENERAL");

  for (i = 0; i < A->n; i++)
    {
      ip = P->pi[i];
      for (j = 0; j < A->m; j++)
	B->mx[i][j] = A->mx[ip][j];
    }
  return eslOK;
}

/* Function:  esl_dmx_LUP_decompose()
 *
 * Purpose:   Calculates a permuted LU decomposition of square matrix
 *            <A>; upon return, <A> is replaced by this decomposition,
 *            where <U> is in the lower triangle (inclusive of the 
 *            diagonal) and <L> is the upper triangle (exclusive of
 *            diagonal, which is 1's by definition), and <P> is the
 *            permutation matrix. Caller provides an allocated 
 *            permutation matrix <P> compatible with the square matrix
 *            <A>.
 *            
 *            Implements Gaussian elimination with pivoting 
 *            \citep[p.~759]{Cormen99}.
 *
 * Throws:    <eslEINVAL> if <A> isn't square, or if <P> isn't the right
 *            size for <A>, or if <A> isn't of general type.
 */
int
esl_dmx_LUP_decompose(ESL_DMATRIX *A, ESL_PERMUTATION *P)
{
  int    i,j,k;
  int    kpiv = 0;    // initialization serves to quiet overzealous static analyzers
  double max;
  double swap;

  if (A->n    != A->m)       ESL_EXCEPTION(eslEINVAL, "matrix isn't square");
  if (P->n    != A->n)       ESL_EXCEPTION(eslEINVAL, "permutation isn't the right size");
  if (A->type != eslGENERAL) ESL_EXCEPTION(eslEINVAL, "matrix isn't of general type");
  esl_permutation_Reuse(P);

  for (k = 0; k < A->n-1; k++)
    {
      /* Identify our pivot; find row with maximum value in col[k]. 
       * This is guaranteed to succeed and set <kpiv> 
       * (no matter what a static analyzer tells you)
       */
      max  = 0.; 
      for (i = k; i < A->n; i++)
	if (fabs(A->mx[i][k]) > max) {
	  max = fabs(A->mx[i][k]);
	  kpiv = i;
	}
      if (max == 0.) ESL_EXCEPTION(eslEDIVZERO, "matrix is singular");
      
      /* Swap those rows (k and kpiv);
       * and keep track of that permutation in P. (misuse j for swapping integers)
       */
      j = P->pi[k]; P->pi[k] = P->pi[kpiv]; P->pi[kpiv] = j;
      for (j = 0; j < A->m; j++)
	{ swap = A->mx[k][j]; A->mx[k][j] = A->mx[kpiv][j]; A->mx[kpiv][j] = swap; }

      /* Gaussian elimination for all rows k+1..n.
       */
      for (i = k+1; i < A->n; i++)
	{
	  A->mx[i][k] /= A->mx[k][k];
	  for (j = k+1; j < A->m; j++)
	    A->mx[i][j] -= A->mx[i][k] * A->mx[k][j];
	}
    }
  return eslOK;
}


/* Function:  esl_dmx_LU_separate()
 *
 * Purpose:   Separate a square <LU> decomposition matrix into its two
 *            triangular matrices <L> and <U>. Caller provides two
 *            allocated <L> and <U> matrices of same size as <LU> for
 *            storing the results.
 *            
 *            <U> may be an upper triangular matrix in either unpacked
 *            (<eslGENERAL>) or packed (<eslUPPER>) form.
 *            <LU> and <L> must be of <eslGENERAL> type.
 *
 * Throws:    <eslEINVAL> if <LU>, <L>, <U> are not of compatible dimensions,
 *            or if <LU> or <L> aren't of general type. 
 */
int
esl_dmx_LU_separate(const ESL_DMATRIX *LU, ESL_DMATRIX *L, ESL_DMATRIX *U)
{
  int i,j;

  if (LU->n    != LU->m)      ESL_EXCEPTION(eslEINVAL, "LU isn't square");
  if (L->n     != L->m)       ESL_EXCEPTION(eslEINVAL, "L isn't square");
  if (U->n     != U->m)       ESL_EXCEPTION(eslEINVAL, "U isn't square");
  if (LU->n    != L->n)       ESL_EXCEPTION(eslEINVAL, "LU, L have incompatible dimensions");
  if (LU->n    != U->n)       ESL_EXCEPTION(eslEINVAL, "LU, U have incompatible dimensions");
  if (LU->type != eslGENERAL) ESL_EXCEPTION(eslEINVAL, "matrix isn't of general type");
  if (L->type  != eslGENERAL) ESL_EXCEPTION(eslEINVAL, "matrix isn't of general type");

  esl_dmatrix_SetZero(L);
  esl_dmatrix_SetZero(U);

  for (i = 0; i < LU->n; i++)
    for (j = i; j < LU->m; j++)
      U->mx[i][j] = LU->mx[i][j];

  for (i = 0; i < LU->n; i++) 
    {
      L->mx[i][i] = 1.;
      for (j = 0; j < i; j++)
	L->mx[i][j] = LU->mx[i][j];
    }
  return eslOK;
}

/* Function:  esl_dmx_Invert()
 *
 * Purpose:   Calculates the inverse of square matrix <A>, and stores the
 *            result in matrix <Ai>. Caller provides an allocated
 *            matrix <Ai> of same dimensions as <A>. Both must be
 *            of type <eslGENERAL>.
 *            
 *            Peforms the inversion by LUP decomposition followed by 
 *            forward/back-substitution \citep[p.~753]{Cormen99}.
 *
 * Throws:    <eslEINVAL> if <A>, <Ai> do not have same dimensions, 
 *                        if <A> isn't square, or if either isn't of
 *                        type <eslGENERAL>.
 *            <eslEMEM>   if internal allocations (for LU, and some other
 *                         bookkeeping) fail.
 */
int
esl_dmx_Invert(const ESL_DMATRIX *A, ESL_DMATRIX *Ai)
{
  ESL_DMATRIX      *LU = NULL;
  ESL_PERMUTATION  *P  = NULL;
  double           *y  = NULL;	/* column vector, intermediate calculation   */
  double           *b  = NULL;	/* column vector of permuted identity matrix */
  int               i,j,k;
  int               status;

  if (A->n     != A->m)                   ESL_EXCEPTION(eslEINVAL, "matrix isn't square");
  if (A->n     != Ai->n || A->m != Ai->m) ESL_EXCEPTION(eslEINVAL, "matrices are different size");
  if (A->type  != eslGENERAL)             ESL_EXCEPTION(eslEINVAL, "matrix A not of general type");
  if (Ai->type != eslGENERAL)             ESL_EXCEPTION(eslEINVAL, "matrix B not of general type");

  /* Copy A to LU, and do an LU decomposition.
   */
  if ((LU = esl_dmatrix_Create(A->n, A->m))    == NULL)  { status = eslEMEM; goto ERROR; }
  if ((P  = esl_permutation_Create(A->n))      == NULL)  { status = eslEMEM; goto ERROR; }
  if (( status = esl_dmatrix_Copy(A, LU))      != eslOK) goto ERROR;
  if (( status = esl_dmx_LUP_decompose(LU, P)) != eslOK) goto ERROR;

  /* Now we have:
   *   PA = LU
   *   
   * to invert a matrix A, we want A A^-1 = I;
   * that's PAx = Pb, for columns x of A^-1 and b of the identity matrix;
   * and that's n equations LUx = Pb;
   * 
   * so, solve Ly = Pb for y by forward substitution;
   * then Ux = y by back substitution;
   * x is then a column of A^-1.
   * 
   * Do that for all columns.
   */
  ESL_ALLOC(b, sizeof(double) * A->n);
  ESL_ALLOC(y, sizeof(double) * A->n);

  for (k = 0; k < A->m; k++)	/* for each column... */
    {
      /* build Pb for column j of the identity matrix */
      for (i = 0; i < A->n; i++)
	if (P->pi[i] == k) b[i] = 1.; else b[i] = 0.;

      /* forward substitution
       */
      for (i = 0; i < A->n; i++)
	{
	  y[i] = b[i];
	  for (j = 0; j < i; j++) y[i] -= LU->mx[i][j] * y[j];
	}

      /* back substitution
       */
      for (i = A->n-1; i >= 0; i--)
	{
	  Ai->mx[i][k] = y[i];
	  for (j = i+1; j < A->n; j++) Ai->mx[i][k] -= LU->mx[i][j] * Ai->mx[j][k];
	  Ai->mx[i][k] /= LU->mx[i][i];
	}
    }

  free(b);
  free(y);
  esl_dmatrix_Destroy(LU);
  esl_permutation_Destroy(P);
  return eslOK;

 ERROR:
  if (y  != NULL) free(y);
  if (b  != NULL) free(b);
  if (LU != NULL) esl_dmatrix_Destroy(LU);
  if (P  != NULL) esl_permutation_Destroy(P);
  return status;
}


/*****************************************************************
 * 7. Optional: interoperability with GSL
 *****************************************************************/
#ifdef HAVE_LIBGSL

#include <gsl/gsl_matrix.h>

int
esl_dmx_MorphGSL(const ESL_DMATRIX *E, gsl_matrix **ret_G)
{
  gsl_matrix *G = NULL;
  int i,j;

  if (E->type != eslGENERAL) ESL_EXCEPTION(eslEINVAL, "can only morph general matrices to GSL right now");

  G = gsl_matrix_alloc(E->m, E->n);
  for (i = 0; i < E->m; i++)
    for (j = 0; j < E->n; j++)
      gsl_matrix_set(G, i, j, E->mx[i][j]);
  *ret_G = G;
  return eslOK;
}

int
esl_dmx_UnmorphGSL(const gsl_matrix *G, ESL_DMATRIX **ret_E)
{
  ESL_DMATRIX *E = NULL;
  int i,j;
  
  if ((E = esl_dmatrix_Create(G->size1, G->size2)) == NULL) return eslEMEM;
  for (i = 0; i < G->size1; i++)
    for (j = 0; j < G->size2; j++)
      E->mx[i][j] = gsl_matrix_get(G, i, j);
  *ret_E = E;
  return eslOK;
}

#endif /*HAVE_LIBGSL*/

/*****************************************************************
 * 8. Optional: Interfaces to LAPACK
 *****************************************************************/
#ifdef HAVE_LIBLAPACK

/* To include LAPACK code, you need to:
 *   1. declare the C interface to the Fortran routine,
 *      appending _ to the Fortran routine's name (dgeev becomes dgeev_)
 *      
 *   2. Remember to transpose matrices into column-major
 *      Fortran form
 *      
 *   3. everything must be passed by reference, not by value
 *   
 *   4. you don't need any include files, just lapack.a
 *   
 *   5. Add -llapack to the compile line.
 *      (It doesn't appear that blas or g2c are needed?)
 */   

/* Declare the C interface to the Fortran77 dgeev routine
 * provided by the LAPACK library:
 */
extern void  dgeev_(char *jobvl, char *jobvr, int *n, double *a,
                    int *lda, double *wr, double *wi, double *vl,
                    int *ldvl, double *vr, int *ldvr,
                    double *work, int *lwork, int *info);


/* Function:  esl_dmx_Diagonalize()
 * Incept:    SRE, Thu Mar 15 09:28:03 2007 [Janelia]
 *
 * Purpose:   Given a square real matrix <A>, diagonalize it:
 *            solve for $U^{-1} A U = diag(\lambda_1... \lambda_n)$.
 *            
 *            Upon return, <ret_Er> and <ret_Ei> are vectors
 *            containing the real and complex parts of the eigenvalues
 *            $\lambda_i$; <ret_UL> is the $U^{-1}$ matrix containing
 *            the left eigenvectors; and <ret_UR> is the $U$ matrix
 *            containing the right eigenvectors.
 *            
 *            <ret_UL> and <ret_UR> are optional; pass <NULL> for
 *            either if you don't want that set of eigenvectors.
 *
 *            This is a C interface to the <dgeev()> routine in the
 *            LAPACK linear algebra library.
 *            
 * Args:      A       -  square nxn matrix to diagonalize
 *            ret_Er  - RETURN: real part of eigenvalues (0..n-1)
 *            ret_Ei  - RETURN: complex part of eigenvalues (0..n-1)
 *            ret_UL  - optRETURN: nxn matrix of left eigenvectors
 *            ret_UR  - optRETURN: 
 *
 * Returns:   <eslOK> on success.
 *            <ret_Er> and <ret_Ei> (and <ret_UL>,<ret_UR> when they are
 *            requested) are allocated here, and must be free'd by the caller.
 *
 * Throws:    <eslEMEM> on allocation failure.
 *            In this case, the four return pointers are returned <NULL>.
 *
 * Xref:      J1/19.
 */
int
esl_dmx_Diagonalize(const ESL_DMATRIX *A, double **ret_Er, double **ret_Ei, 
		    ESL_DMATRIX **ret_UL, ESL_DMATRIX **ret_UR)
{
  int          status;
  double      *Er   = NULL;
  double      *Ei   = NULL;
  ESL_DMATRIX *At   = NULL;
  ESL_DMATRIX *UL   = NULL;
  ESL_DMATRIX *UR   = NULL;
  double      *work = NULL;
  char   jobul, jobur;
  int    lda;
  int    ldul, ldur;
  int    lwork;
  int    info;

  if (A->n != A->m) ESL_EXCEPTION(eslEINVAL, "matrix isn't square");

  if ((At = esl_dmatrix_Clone(A))          == NULL)       { status = eslEMEM; goto ERROR; } 
  if ((UL = esl_dmatrix_Create(A->n,A->n)) == NULL)       { status = eslEMEM; goto ERROR; }
  if ((UR = esl_dmatrix_Create(A->n,A->n)) == NULL)       { status = eslEMEM; goto ERROR; }
  ESL_ALLOC(Er,   sizeof(double) * A->n);
  ESL_ALLOC(Ei,   sizeof(double) * A->n);
  ESL_ALLOC(work, sizeof(double) * 8 * A->n);

  jobul = (ret_UL == NULL) ? 'N' : 'V';	/* do we want left eigenvectors? */
  jobur = (ret_UR == NULL) ? 'N' : 'V'; /* do we want right eigenvectors? */
  lda   = A->n; 
  ldul  = A->n;
  ldur  = A->n;
  lwork = 8*A->n;

  /* Fortran convention is colxrow, not rowxcol; so transpose
   * a copy of A before passing it to a Fortran routine.
   */
  esl_dmx_Transpose(At);

  /* The Fortran77 interface call to LAPACK's dgeev().
   * All args must be passed by reference.
   * Fortran 2D arrays are 1D: so pass the A[0] part of a DSMX.
   */
  dgeev_(&jobul, &jobur, &(At->n), At->mx[0], &lda, Er, Ei, 
	 UL->mx[0], &ldul, UR->mx[0], &ldur, work, &lwork, &info);
  if (info < 0) ESL_XEXCEPTION(eslEINVAL, "argument %d to LAPACK dgeev is invalid", -info);
  if (info > 0) ESL_XEXCEPTION(eslEINVAL, 
			       "diagonalization failed; only eigenvalues %d..%d were computed",
			       info+1, At->n);

  /* Now, UL, UR are transposed (col x row), so transpose them back to
   * C language convention.
   */
  esl_dmx_Transpose(UL);
  esl_dmx_Transpose(UR);

  esl_dmatrix_Destroy(At);
  if (ret_UL != NULL) *ret_UL = UL; else esl_dmatrix_Destroy(UL);
  if (ret_UR != NULL) *ret_UR = UR; else esl_dmatrix_Destroy(UR);
  if (ret_Er != NULL) *ret_Er = Er; else free(Er);
  if (ret_Ei != NULL) *ret_Ei = Ei; else free(Ei);
  free(work);
  return eslOK;

 ERROR:
  if (ret_UL != NULL) *ret_UL = NULL;
  if (ret_UR != NULL) *ret_UR = NULL;
  if (ret_Er != NULL) *ret_Er = NULL;
  if (ret_Ei != NULL) *ret_Ei = NULL;
  if (At   != NULL) esl_dmatrix_Destroy(At);
  if (UL   != NULL) esl_dmatrix_Destroy(UL);
  if (UR   != NULL) esl_dmatrix_Destroy(UR);
  if (Er   != NULL) free(Er);
  if (Ei   != NULL) free(Ei);
  if (work != NULL) free(work);
  return status;
}


#endif /*HAVE_LIBLAPACK*/

/*****************************************************************
 * 9. Unit tests
 *****************************************************************/ 
#ifdef eslDMATRIX_TESTDRIVE

static void 
utest_misc_ops(void)
{
  char *msg = "miscellaneous unit test failed";
  ESL_DMATRIX *A, *B, *C;
  int  n = 20;

  if ((A = esl_dmatrix_Create(n,n)) == NULL) esl_fatal(msg);
  if ((B = esl_dmatrix_Create(n,n)) == NULL) esl_fatal(msg);
  if ((C = esl_dmatrix_Create(n,n)) == NULL) esl_fatal(msg);
  
  if (esl_dmatrix_SetIdentity(A)    != eslOK) esl_fatal(msg);   /* A=I */
  if (esl_dmx_Invert(A, B)          != eslOK) esl_fatal(msg);	/* B=I^-1=I */
  if (esl_dmx_Multiply(A,B,C)       != eslOK) esl_fatal(msg);	/* C=I */
  if (esl_dmx_Transpose(A)          != eslOK) esl_fatal(msg);   /* A=I still */

  if (esl_dmx_Scale(A, 0.5)         != eslOK) esl_fatal(msg);	/* A= 0.5I */
  if (esl_dmx_AddScale(B, -0.5, C)  != eslOK) esl_fatal(msg);	/* B= 0.5I */
  
  if (esl_dmx_Add(A,B)              != eslOK) esl_fatal(msg);	/* A=I */
  if (esl_dmx_Scale(B, 2.0)         != eslOK) esl_fatal(msg);	/* B=I */
  
  if (esl_dmatrix_Compare(A, B, 1e-12) != eslOK) esl_fatal(msg);
  if (esl_dmatrix_Compare(A, C, 1e-12) != eslOK) esl_fatal(msg);
  if (esl_dmatrix_Copy(B, C)           != eslOK) esl_fatal(msg);
  if (esl_dmatrix_Compare(A, C, 1e-12) != eslOK) esl_fatal(msg);

  esl_dmatrix_Destroy(A);
  esl_dmatrix_Destroy(B);
  esl_dmatrix_Destroy(C);
  return;
}


static void
utest_Invert(ESL_DMATRIX *A)
{
  char *msg = "Failure in matrix inversion unit test";
  ESL_DMATRIX *Ai = NULL;
  ESL_DMATRIX *B  = NULL;
  ESL_DMATRIX *I  = NULL;

  if ((Ai = esl_dmatrix_Create(A->n, A->m)) == NULL)  esl_fatal(msg);
  if ((B  = esl_dmatrix_Create(A->n, A->m)) == NULL)  esl_fatal(msg);
  if ((I  = esl_dmatrix_Create(A->n, A->m)) == NULL)  esl_fatal(msg);
  if (esl_dmx_Invert(A, Ai)                 != eslOK) esl_fatal("matrix inversion failed");
  if (esl_dmx_Multiply(A,Ai,B)              != eslOK) esl_fatal(msg);
  if (esl_dmatrix_SetIdentity(I)            != eslOK) esl_fatal(msg);
  if (esl_dmatrix_Compare(B,I, 1e-12)       != eslOK) esl_fatal("inverted matrix isn't right");
  
  esl_dmatrix_Destroy(Ai);
  esl_dmatrix_Destroy(B);
  esl_dmatrix_Destroy(I);
  return;
}


#endif /*eslDMATRIX_TESTDRIVE*/



/*****************************************************************
 * 10. Test driver
 *****************************************************************/ 

/*   gcc -g -Wall -o test -I. -L. -DeslDMATRIX_TESTDRIVE esl_dmatrix.c -leasel -lm
 */
#ifdef eslDMATRIX_TESTDRIVE
#include "easel.h"
#include "esl_dmatrix.h"
#include "esl_random.h"

int main(void)
{
  ESL_RANDOMNESS *r;
  ESL_DMATRIX *A;
  int          n    = 30;
  int          seed = 42;
  int          i,j;
  double       range = 100.;

  /* Create a square matrix with random values from  -range..range */
  if ((r = esl_randomness_Create(seed)) == NULL) esl_fatal("failed to create random source");
  if ((A = esl_dmatrix_Create(n, n))    == NULL) esl_fatal("failed to create matrix");
  for (i = 0; i < n; i++)
    for (j = 0; j < n; j++)
      A->mx[i][j] = esl_random(r) * range * 2.0 - range;

  utest_misc_ops();
  utest_Invert(A);

  esl_randomness_Destroy(r);
  esl_dmatrix_Destroy(A);
  return 0;
}
#endif /*eslDMATRIX_TESTDRIVE*/


/*****************************************************************
 * 11. Examples
 *****************************************************************/ 

/*   gcc -g -Wall -o example -I. -DeslDMATRIX_EXAMPLE esl_dmatrix.c easel.c -lm
 */
#ifdef eslDMATRIX_EXAMPLE
/*::cexcerpt::dmatrix_example::begin::*/
#include "easel.h"
#include "esl_dmatrix.h"

int main(void)
{
  ESL_DMATRIX *A, *B, *C;

  A = esl_dmatrix_Create(4,4);
  B = esl_dmatrix_Create(4,4);
  C = esl_dmatrix_Create(4,4);
  
  esl_dmatrix_SetIdentity(A);
  esl_dmatrix_Copy(A, B);

  esl_dmx_Multiply(A,B,C);

  esl_dmatrix_Dump(stdout, C, NULL, NULL);

  esl_dmatrix_Destroy(A);
  esl_dmatrix_Destroy(B);
  esl_dmatrix_Destroy(C);
  return 0;
}
/*::cexcerpt::dmatrix_example::end::*/
#endif /*eslDMATRIX_EXAMPLE*/