/* See SOMatrix.h */
/* Last edited on 2003-01-15 21:11:59 by anamaria */

#include <SOMatrix.h>

#include <SOBasic.h>
#include <vec.h>
#include <filefmt.h>
#include <fget.h>
#include <nget.h>
#include <js.h>
#include <stdio.h>
#include <stdlib.h>
#include <math.h>
#include <limits.h>
  
void SOMatrix_UnpackRow(MatEntry_vec *Aents, nat row, double_vec Ai, nat *kA);
  /* Unpacks row {i} of the matrix {A} as an ordinary vector of double
    {Ai}. On input, the variable {*kA} must be the index into {Aents ==}
{    A.ents} of the first element of that row. On output, {*kA} will be
    the index of the first element of row {i+1}. */

void SOMatrix_UnpackHalfRow(MatEntry_vec *Lents, nat row, double_vec Li, nat *kL);
  /* Same as {UnpackRow}, but expects and checks that {L} is lower triangular ---
    meaning that there are no elements {L[i,j]} with {j > i}. */

void SOMatrix_PackRow(double_vec Ci, nat row, MatEntry_vec *Cents, nat *nC);
  /* Appends the non-zero entries of {Ci} as one more row of elements
    of a matrix {C}. Will reallocate the element list {Cents == C.ents} if
    necessary. On input, the variable {*nC} must be the count of
    elements already stored in {Cents^}; this count will be updated on
    exit. Don't forget to trim the element list after adding the last
    row. */

SOMatrix SOMatrix_Null(nat rows, nat cols)
  {     
    MatEntry_vec ents = MatEntry_vec_new(0);
    return (SOMatrix){rows, rows, ents};
  }

SOMatrix SOMatrix_BinaryOp
  ( SOMatrix A, SOMatrix B,
    double func(int row, int col, double va, double vb),
    bool skip00
  )
  { MatEntry *Ael = A.ents.el; int nA = A.ents.nel;
    MatEntry *Bel = B.ents.el; int nB = B.ents.nel;
    MatEntry_vec Cents = MatEntry_vec_new(nA + nB);
    nat kA = 0, kB = 0, nC = 0;
    nat row, col;
    double va, vb, vc;
    assert(A.rows == B.rows, "inconsistent rows");
    assert(A.cols == B.cols, "inconsistent cols");
    row = -1; col = A.cols-1;
    while (TRUE)
      { /* Determine indices {[row,col]} of next element to be computed: */
        if (! skip00)
          { /* Must evaluate even elements which are missing in {A} and {B}: */
            col++; 
            if (col >= A.cols) { row++; col = 0; }
          }
        else if ((kA >= nA) && (kB >= nB))
          { row = A.rows; col = 0; }
        else if (kA >= nA)
          { row = Bel[kB].row; col = Bel[kB].col; }
        else if (kB >= nB)
          { row = Ael[kA].row; col = Ael[kA].col; }
        else if (Ael[kA].row == Bel[kB].row)
          { if (Ael[kA].col <= Bel[kB].col)
              { row = Ael[kA].row; col = Ael[kA].col; }
            else
              { row = Bel[kB].row; col = Bel[kB].col; }
          }
        else if (Ael[kA].row < Bel[kB].row)
          { row = Ael[kA].row; col = Ael[kA].col; }
        else
          { row = Bel[kB].row; col = Bel[kB].col; }
        
        /* Are we done already? */
        if (row >= A.rows) { break; }
        
        /* Get element values from {A} and {B}, and skip over them: */
        if ((kA < nA) && (row == Ael[kA].row) && (col == Ael[kA].col))
          { va = Ael[kA].va; kA++; }
        else
          { va = 0.0; }
        if ((kB < nB) && (row == Bel[kB].row) && (col == Bel[kB].col))
          { vb = Bel[kB].va; kB++; }
        else
          { vb = 0.0; }
        
        /* Compute element of result: */
        if (skip00 && (va == 0.0) && (vb == 0.0))
          { vc = 0.0; } 
        else
          { vc = func(row, col, va, vb); }

        /* Append element to {C}: */
        if (vc != 0.0)
          { MatEntry_vec_expand(&Cents, nC);
            { MatEntry *Cij = &(Cents.el[nC]);
              Cij->row = row; Cij->col = col;
              Cij->va = vc; nC++;
            }
          }
      }
    MatEntry_vec_trim(&Cents, nC);
    return (SOMatrix){A.rows, B.rows, Cents};
  }

SOMatrix SOMatrix_Mix(double alpha, SOMatrix A, double beta, SOMatrix B)
  { auto double mix(int row, int col, double va, double vb);
    
    double mix(int row, int col, double va, double vb)
      { return alpha*va + beta*vb; }
      
    return SOMatrix_BinaryOp(A, B, mix, TRUE);
  }
  
SOMatrix SOMatrix_RelDiff(SOMatrix A, SOMatrix B)
  { auto double rdiff(int row, int col, double va, double vb);
    
    double rdiff(int row, int col, double va, double vb)
      { double ma = fabs(va), mb = fabs(vb);
        double m = (ma > mb ? ma : mb);
        if (m == 0.0) 
          { return 0.0; }
        else
          { double ra = va/m, rb = vb/m;
            return 0.5*(rb - ra)/sqrt(ra*ra + rb*rb);
          }
      }
      
    return SOMatrix_BinaryOp(A, B, rdiff, TRUE);
  }

SOMatrix SOMatrix_Scale(double alpha, SOMatrix A)
  { MatEntry *Ael = A.ents.el; int nA = A.ents.nel;
    nat nC = 0, kA;
    MatEntry_vec Cents = MatEntry_vec_new((alpha == 0.0 ? 0 : nA));
    for (kA = 0; kA < nA; kA++)
      { MatEntry *Aij = &(Ael[kA]);
        MatEntry *Cij = &(Cents.el[nC]);
        double va = alpha * Aij->va;
        if (va != 0.0) { (*Cij) = (*Aij); Cij->va = va; nC++; }
      }
    MatEntry_vec_trim(&Cents, nC);
    return (SOMatrix){A.rows, A.cols, Cents};
  }

void SOMatrix_GetDiagonal(SOMatrix A, double_vec d)
  { MatEntry *Ael = A.ents.el; int nA = A.ents.nel;
    int nd = (A.rows < A.cols ? A.rows : A.cols);
    int kd = 0, kA;
    assert(d.nel == nd, "d has wrong size");
    for (kA = 0; kA < nA; kA++)
      { MatEntry *Aij = &(Ael[kA]);
        if (Aij->row == Aij->col) 
          { int i = Aij->row;
            assert(i <= kd, "matrix entries out of order");
            while (kd < i) { d.el[kd] = 0.0; kd++; }
            d.el[kd] = Aij->va;
            kd++;
          }
      }
    while (kd < nd) { d.el[kd] = 0.0; kd++; }
  }

SOMatrix SOMatrix_Mul(SOMatrix A, SOMatrix B)
  { int nA = A.ents.nel;
    int nB = B.ents.nel;
    nat maxnel = (nA > nB ? nA : nB);
    MatEntry_vec Cents = MatEntry_vec_new(maxnel);
    double_vec a = double_vec_new(A.rows);
    double_vec c = double_vec_new(B.cols);
    int i; nat kA = 0, nC = 0;

    assert(A.cols == B.rows , "incompatible sizes");
    for (i = 0; i < A.rows; i++)
      { SOMatrix_UnpackRow(&A.ents, i, a, &kA);
        SOMatrix_MulRow(a, B, c);
        SOMatrix_PackRow(c, i, &Cents, &nC);
      }
    MatEntry_vec_trim(&Cents, nC);
    free(a.el); free(c.el);
    return (SOMatrix){A.rows, B.cols, Cents};
  }

void SOMatrix_MulRow(double_vec x, SOMatrix A, double_vec y)
  { MatEntry *Ael = A.ents.el; int nA = A.ents.nel;
    int j, kA;
    assert(x.nel == A.rows, "vector has wrong size");
    assert(y.nel == A.cols, "vector has wrong size");
    for (j = 0; j < A.cols; j++) { y.el[j] = 0.0; }
    for (kA = 0; kA < nA; kA++)
      { MatEntry *Aij = &(Ael[kA]);
        nat i = Aij->row, j = Aij->col;
        y.el[j] += x.el[i] * Aij->va;
      }
  }

void SOMatrix_MulCol(SOMatrix A, double_vec x, double_vec y)
  { MatEntry *Ael = A.ents.el; int nA = A.ents.nel;
    int i, kA;
    assert(x.nel == A.cols, "vector has wrong size");
    assert(y.nel == A.rows, "vector has wrong size");
    for (i = 0; i < A.cols; i++) { y.el[i] = 0.0; }
    for (kA = 0; kA < nA; kA++)
      { MatEntry *Aij = &(Ael[kA]);
        nat i = Aij->row, j = Aij->col;
        y.el[i] += Aij->va * x.el[j];
      }
  }

void SOMatrix_DivRow(double_vec y, SOMatrix L, double_vec x)
  { MatEntry *Lel = L.ents.el; int nL = L.ents.nel;
    int i, kL;

    assert(L.rows == L.cols, "matrix is not square");
    assert(y.nel == L.cols, "vector has wrong size");
    assert(x.nel == L.cols, "vector has wrong size");

    for (i = 0; i < L.cols; i++) { x.el[i] = y.el[i]; }
    for (kL = nL-1; kL >= 0; kL--)
      { MatEntry *Lij = &(Lel[kL]);
          nat i = Lij->row, j = Lij->col;
          assert(i >= j, "not lower triangular");
          if (i > j)
            { x.el[j] -= Lij->va * x.el[i]; }
          else
            { x.el[j] /= Lij->va; }
        }
    }

void SOMatrix_DivCol(SOMatrix L, double_vec y, double_vec x)
  { int i, j, kL;
    double_vec Li = double_vec_new(L.cols);
    
    assert(L.rows == L.cols, "matrix is not square");
    assert(y.nel == L.cols, "vector has wrong size");
    assert(x.nel == L.cols, "vector has wrong size");
    
    kL = 0;
    for (i = 0; i < L.rows; i++)
      { SOMatrix_UnpackHalfRow(&L.ents, i, Li, &kL);
        x.el[i] = y.el[i];
        for (j = 0; j < i; j++)
          { x.el[i] -= Li.el[j]*x.el[j]; }
        x.el[i] /= Li.el[i];
      }
    free(Li.el);
  }

SOMatrix SOMatrix_PosDiv(SOMatrix L, SOMatrix A)
  { int i, kA, nX;
    MatEntry_vec Xents = MatEntry_vec_new(A.ents.nel);
    double_vec Ai = double_vec_new(A.cols);
    double_vec Xi = double_vec_new(A.cols);
    
    assert(L.rows == L.cols, "L matrix is not square");
    assert(A.cols == L.rows, "incompatible matrix sizes");
    
    kA = 0; nX = 0;
    for (i = 0; i < A.rows; i++)
      { SOMatrix_UnpackRow(&(A.ents), i, Ai, &kA);
        SOMatrix_DivRow(Ai, L, Xi);
        SOMatrix_PackRow(Xi, i, &Xents, &nX);
      }
    free(Ai.el); free(Xi.el);
    MatEntry_vec_trim(&Xents, nX);
    return (SOMatrix){A.rows, A.cols, Xents}; 
  }

SOMatrix SOMatrix_PreDiv(SOMatrix L, SOMatrix A)
  { int i, kL, kA, kX, nX;
    MatEntry_vec Xents = MatEntry_vec_new(A.ents.nel);
    double_vec Li = double_vec_new(L.cols);
    double_vec Xi = double_vec_new(A.cols);
    
    assert(L.rows == L.cols, "L matrix is not square");
    assert(A.rows == L.cols, "A matrix has wrong rows");
    
    kL = 0; kA = 0; nX = 0;
    for (i = 0; i < L.rows; i++)
      { SOMatrix_UnpackRow(&(L.ents), i, Li, &kL);
        SOMatrix_UnpackRow(&(A.ents), i, Xi, &kA);
        kX = 0;
        while(kX < nX)
          { MatEntry *Xjk = &(Xents.el[kX]);
            int j = Xjk->row, k = Xjk->col;
            Xi.el[k] -= Li.el[j] * Xjk->va;
          }
        { int k; for (k = 0; k < A.cols; k++) { Xi.el[k] /= Li.el[i]; } }
        SOMatrix_PackRow(Xi, i, &Xents, &nX);
      }
    free(Li.el); free(Xi.el);
    MatEntry_vec_trim(&Xents, nX);
    return (SOMatrix){A.rows, A.cols, Xents}; 
  }

SOMatrix SOMatrix_Cholesky(SOMatrix A, double minVal)
  { 
    /* Andre-Louis Cholesky (spelled with a 'y'), born in France in
      1875, was a geodesist in the French military. He developed his
      computational procedure to solve geodetic problems, among which
      was, in the words of his eulogy, the /problem of levelling in
      Morocco./ He died in battle in 1918. Benoit published the
      computational method in /Bulletin geodesique/ in 1924.

      The most likely French pronunciation of Cholesky is that heard
      in 'Chopin' -- the English 'Sh' sound. 
        -- David Pattison, Washington, DC */
  
    MatEntry_vec Lents = MatEntry_vec_new(A.rows);
    int i, j, nL = 0, kL, kA = 0;
    double Lii;

    double_vec Ai = double_vec_new(A.rows);
    double_vec Li = double_vec_new(A.rows);
    double_vec dg = double_vec_new(A.rows); /* Diagonal elements */
    
    assert(A.rows == A.cols , "matrix must be square");
    for (i = 0; i < A.rows; i++)
      { SOMatrix_UnpackRow(&A.ents, i, Ai, &kA);
        /* Compute {L[i,j] = X[i,j] - SUM{ L[i,k]*A[j,k] : k = 0..j-1 })/L[j,j]} */
        /* where {X[i,j] = (j <= i ? A[i,j] : 0)} */
        kL = 0;
        for (j = 0; j < i; j++)
          { double sum = 0.0, corr = 0.0;
            while (kL < nL)
              { MatEntry *Lrs = &(Lents.el[kL]); 
                int k = Lrs->col;
                assert(Lrs->row == j, "matrix indices buggy");
                kL++;
                if (k == j)
                  { double Ljj = Lrs->va;
                    assert(Ljj != 0.0, "zero element in diagonal");
                    Li.el[j] = (Ai.el[j] - sum)/Ljj;
                    break;
                  } 
                else
                  { double Ljk = Lrs->va;
                    double term = Li.el[k]*Ljk;
                    /* Kahan's summation formula: */
                    double tcorr = term - corr;
                    double newSum = sum + tcorr;
                    corr = (newSum - sum) - tcorr;
                    sum = newSum;
                    assert (k < j, "Cholesky's upper half nonzero");
                  }
              }
          }
        
        /* Compute {L[i,i] = sqrt(L[i,i] - SUM{ L[i,j]^2 : j = 0..i-1 })} */
        { double sum = 0.0, corr = 0.0;
          for (j = 0; j < i; j++)
            { double w = Li.el[j], term = w*w; 
              /* Kahan's summation formula: */
              double tcorr = term - corr;
              double newSum = sum + tcorr;
              corr = (newSum - sum) - tcorr;
              sum = newSum;
            }
          Lii = Ai.el[i] - sum;
        }
        if (Lii < 0.0)
          { fprintf(stderr, "matrix is not positive definite?\n");
            fprintf(stderr, "row %04d Lii = %24.15e\n", i, Lii);
            for (j = 0; j <= i; j++)
              { if (Li.el[j] != 0.0) 
                  { fprintf(stderr, "  %04d %24.15e\n", j, Li.el[j]); }
              }
            fprintf(stderr, "\n");
            assert(FALSE, "aborted");
          }
        Li.el[i] = sqrt(Lii);
        /* Save diagonal elements: */
        dg.el[i] = Li.el[i];
        /* Cleanup small entries: */
        for (j = 0; j < i; j++)
          { if ((minVal > 0.0) && (fabs(Li.el[j]) < minVal*sqrt(dg.el[i]*dg.el[j])))
              { Li.el[j] = 0.0; }
          }
        for (j = i+1; j < A.cols; j++) { Li.el[j] = 0.0; }
        SOMatrix_PackRow(Li, i, &Lents, &nL);
      }
    free(Ai.el); free(Li.el); free(dg.el);
    MatEntry_vec_trim(&Lents, nL);

    return (SOMatrix){A.rows, A.cols, Lents};
  }
  
void SOMatrix_GaussSeidel(SOMatrix A, double_vec b, double omega, double_vec x)
  { MatEntry *Ael = A.ents.el; int nA = A.ents.nel;
    
    int i, kA = 0;
    double Aii;

    assert(A.rows == A.cols, "matrix must be square");
    assert(b.nel == A.cols, "vector has wrong size");
    assert(x.nel == A.cols, "vector has wrong size");
    for (i = 0; i < A.rows; i++)
      { double s = 0.0;
        while ((kA < nA) && (Ael[kA].row  == i))
          { MatEntry *Aij = &(Ael[kA]);
            nat j = Aij->col;
            if (j == i)
              { Aii = Aij->va; }
            else
              { s += Aij->va * x.el[j]; }
            kA++;
          }
        x.el[i] = (1.0 - omega) * x.el[i] + omega*(b.el[i] - s)/Aii;
      }
  }

double SOMatrix_ConjugateGradient(SOMatrix A, double_vec b, double_vec x)
  {
    double rsq;
    double_vec xi = double_vec_new(A.cols);
    double_vec xj = double_vec_new(A.cols);
    double_vec g = double_vec_new(A.cols);
    double_vec h = double_vec_new(A.cols);
    double eps = 1.0e-6;

    auto bool ConjugateGradientIteration(void);
      /* Performs one `macro' iteration of the Conjugate Gradient
        algorithm for solving the linear system {A x == b}.
        Returns TRUE iff the iteration converged. */  
    
    bool ConjugateGradientIteration(void) 
      {
        double rp = 0.0, bsq = 0.0;
        double anum, aden, gam, dgg, gg;
        int iter, j;
        SOMatrix_MulCol(A, x, xi);
        for (j = 0; j < A.cols; j++)
          { bsq += b.el[j]*b.el[j];
            xi.el[j] -= b.el[j];
            rp += xi.el[j]*xi.el[j];
          }
        SOMatrix_MulRow(xi, A, g);
        for (j = 0; j < A.cols; j++)
          { g.el[j] = -g.el[j];
            h.el[j] = g.el[j];
          }
        for (iter = 1; iter <= 10*A.cols; iter++)
          { fprintf(stderr, "*");
            SOMatrix_MulCol(A, h, xi);
            anum = 0.0;
            aden = 0.0;
            for (j = 0; j < A.cols; j++)
              { anum = anum + g.el[j]*h.el[j];
                aden = aden + xi.el[j]*xi.el[j];
              }
            assert(aden != 0.0, "divide by zero"); 
            anum = anum/aden;
            for (j = 0; j < A.cols; j++)
              { xi.el[j] = x.el[j];
                x.el[j] += anum * h.el[j];
              }
            SOMatrix_MulCol(A, x, xj);
            rsq = 0.0;
            for (j = 0; j < A.cols; j++)
              { xj.el[j] -= b.el[j];
                rsq += xj.el[j]*xj.el[j];
              }
            if ((rsq == rp) || (rsq <= bsq*((double)A.cols)*eps*eps))
              { fprintf(stderr, "1");
                return TRUE ;
              }
            if (rsq > rp)
              { x = xi;
                fprintf(stderr, "2");
                return FALSE ;
              }
            rp = rsq;
            SOMatrix_MulRow(xj, A, xi);
            gg = 0.0;
            dgg = 0.0;
            for (j = 0; j < A.cols; j++)
              { gg = gg + g.el[j]*g.el[j];
                dgg = dgg + (xi.el[j] + g.el[j])*xi.el[j];
              }
            if (gg == 0.0)
              { fprintf(stderr, "3");
                return TRUE;
              }
            gam = dgg/gg;
            for (j = 0; j < A.cols; j++)
              { g.el[j] = -xi.el[j];
                h.el[j] = g.el[j] + gam*h.el[j];
              }
          }
        fprintf(stderr, "4");
        return TRUE;
      }

    int i;
    
    assert(A.rows == A.cols , "matrix must be square");
    assert(b.nel == A.cols, "wrong dimensions");
    assert(x.nel == A.cols, "wrong dimensions");
    for (i = 1; i <= 3; i++)
      { fprintf(stderr, "!");
        if (ConjugateGradientIteration()) { return rsq; }
      }
    free(xi.el); free(xj.el);
    free(g.el); free(h.el);
    return rsq;
  }

nat SOMatrix_Find(MatEntry *r, nat n, nat row, nat col)
  { /* Binary search in the element list. */
    int p = 0, q = n-1;
    while (p <= q)
      { int m = (p + q)/2;
        MatEntry *rm = &(r[m]);
        if ((rm->row < row) || ((rm->row == row) && (rm->col < col)))
          { p = m + 1; }
        else if ((rm->row > row) || ((rm->row == row) && (rm->col > col)))
          { q = m - 1; }
        else
          { return m; }
      }
    return INT_MAX;
  }

double SOMatrix_GetValue(SOMatrix A, nat row, nat col)
  {
    nat k = SOMatrix_Find(A.ents.el, A.ents.nel, row, col);
    return (k == INT_MAX ? 0.0 : A.ents.el[k].va);
  }

#define SOMatrix_FileFormat "2002-11-20"
#define SOListMatrix_FileFormat "2004-09-07"
    
void SOMatrix_Write(FILE *wr, SOMatrix A)
  { MatEntry *Ael = A.ents.el; int nA = A.ents.nel;
    int kA;
    filefmt_write_header(wr, "SOMatrix", SOMatrix_FileFormat);
    fprintf(wr, "rows = %d\n", A.rows);
    fprintf(wr, "cols = %d\n", A.cols);
    fprintf(wr, "elems = %d\n", nA);
    for (kA = 0; kA < nA; kA++) 
      { MatEntry *Aij = &(Ael[kA]);
        fprintf(wr, "%d %d %24.16e\n", Aij->row, Aij->col, Aij->va);
      }
    filefmt_write_footer(wr, "SOMatrix");
    fflush(wr);
  }

SOMatrix SOMatrix_Read(FILE *rd)
  { SOMatrix A;
    int nA;
    filefmt_read_header(rd, "SOMatrix", SOMatrix_FileFormat);
    A.rows = nget_int(rd, "rows"); fget_eol(rd); 
    A.cols = nget_int(rd, "cols"); fget_eol(rd);
    nA = nget_int(rd, "elems"); fget_eol(rd);
    { MatEntry_vec Aents = MatEntry_vec_new(nA);
      int kA;
      for (kA = 0; kA < nA; kA++)
        { Aents.el[kA].row = fget_int(rd); 
          Aents.el[kA].col = fget_int(rd); 
          Aents.el[kA].va = fget_double(rd);
          fget_eol(rd);
        }
      A.ents = Aents;
    }  
    filefmt_read_footer(rd, "SOMatrix");
    return A;
  }

void SOListMatrix_Write(FILE *wr, SOListMatrix A)
  { ListMatEntry *Ael = A.ents.el; int nA = A.ents.nel;
    int kA, lA;
    filefmt_write_header(wr, "SOListMatrix", SOListMatrix_FileFormat);
    fprintf(wr, "rows = %d\n", A.rows);
    fprintf(wr, "cols = %d\n", A.cols);
    fprintf(wr, "length = %d\n", A.length);
    fprintf(wr, "elems = %d\n", nA);
    for (kA = 0; kA < nA; kA++) 
      { ListMatEntry *Aij = &(Ael[kA]);
        fprintf(wr, "%d %d ", Aij->row, Aij->col);

        for (lA = 0; lA < A.length; lA++) 
          fprintf(wr, "%24.16e ", Aij->va[lA]);

        fprintf(wr, "\n");
      }
    filefmt_write_footer(wr, "SOListMatrix");
    fflush(wr);
  }

SOListMatrix SOListMatrix_Read(FILE *rd)
  { SOListMatrix A;
    int nA;
    filefmt_read_header(rd, "SOListMatrix", SOListMatrix_FileFormat);
    A.rows = nget_int(rd, "rows"); fget_eol(rd); 
    A.cols = nget_int(rd, "cols"); fget_eol(rd);
    A.length = nget_int(rd, "length"); fget_eol(rd);
    nA = nget_int(rd, "elems"); fget_eol(rd);
    { ListMatEntry_vec Aents = ListMatEntry_vec_new(nA);
      int kA, lA;
      for (kA = 0; kA < nA; kA++)
        { Aents.el[kA].row = fget_int(rd); 
          Aents.el[kA].col = fget_int(rd); 

          Aents.el[kA].va = (double *)malloc(sizeof(double)*A.length);
          for (lA = 0; lA < A.length; lA++)
            Aents.el[kA].va[lA] = fget_double(rd);

          fget_eol(rd);
        }
      A.ents = Aents;
    }  
    filefmt_read_footer(rd, "SOListMatrix");
    return A;
  }

void SOMatrix_UnpackRow(MatEntry_vec *Aents, nat row, double_vec Ai, nat *kA)
  { int j;
    MatEntry *Aij;
    for (j = 0; j < Ai.nel; j++) { Ai.el[j] = 0.0; }
    while (((*kA) < Aents->nel) && ((Aij = &(Aents->el[(*kA)]))->row == row))
      { Ai.el[Aij->col] = Aij->va; (*kA)++; }
  }

void SOMatrix_UnpackHalfRow(MatEntry_vec *Lents, nat row, double_vec Li, nat *kL)
  { int j;
    MatEntry *Lij;
    for (j = 0; j < Li.nel; j++) { Li.el[j] = 0.0; }
    while (((*kL) < Lents->nel) && ((Lij = &(Lents->el[(*kL)]))->row == row))
      { Li.el[Lij->col] = Lij->va;
        assert(Lij->col <= row, "matrix is not lower triangular");
        (*kL)++;
      }
  }

void SOMatrix_PackRow(double_vec Ci, nat row, MatEntry_vec *Cents, nat *nC)
  { int j;
    for (j = 0; j < Ci.nel; j++)
      { if (Ci.el[j] != 0.0)
          { MatEntry_vec_expand(Cents, *nC);
            { MatEntry *Cij = &(Cents->el[*nC]);
              Cij->row = row; Cij->col = j;
              Cij->va = Ci.el[j];
              (*nC)++;
            }
          }
      }
  }

/* Arrays of {MatEntry}: */

MatEntry_vec MatEntry_vec_new(nat nel)
  { /* This is not a macro only because gcc does not allow cast of struct: */
    vec_t v = vec_new(nel, sizeof(MatEntry));
    MatEntry_vec r;
    r.nel = v.nel; r.el = (MatEntry*)v.el;
    return r;
  }

/* Arrays of {ListMatEntry}: */

ListMatEntry_vec ListMatEntry_vec_new(nat nel)
  { /* This is not a macro only because gcc does not allow cast of struct: */
    vec_t v = vec_new(nel, sizeof(ListMatEntry));
    ListMatEntry_vec r;
    r.nel = v.nel; r.el = (ListMatEntry*)v.el;
    return r;
  }