/* See SOBasisMatrices.h */
/* Last edited on 2003-04-28 19:21:18 by ra014852 */

#include <SOBasisMatrices.h>
#include <SOFunction.h>
#include <SOMatrix.h>
#include <SOBasic.h>
#include <vec.h>
#include <math.h>
#include <limits.h>
#include <js.h>
 
SOMatrix SOBasisMatrices_BuildMatrixGen
  ( Basis F,
    Basis G, 
    void dot(SOFunction *f, SOFunction *g, double *iv), //double dot(SOFunction *f, SOFunction *g), 
    double minVal,
    bool symmetric,
    bool verbose
  )
  { double fg;
    nat m = F.nel;
    nat n = G.nel;
    MatEntry_vec Ments = MatEntry_vec_new(F.nel+G.nel);
    double_vec fsqr = double_vec_new(F.nel);
    double_vec gsqr = double_vec_new(G.nel);
    int i, j; nat kM = 0;
    bool veryverbose = FALSE;
    
    if (minVal > 0.0) 
      { /* Compute the element norms squared, for small entry cleanup: */
        if (verbose) { fprintf(stderr, "computing norms...\n"); }
        for (i = 0; i < m; i++) 
          { dot(F.el[i], F.el[i], &(fsqr.el[i])); } //{ fsqr.el[i] = dot(F.el[i], F.el[i]); }
        for (j = 0; j < n; j++) 
          { //{ gsqr.el[j]=(symmetric?fsqr.el[j]:dot(G.el[j],G.el[j]));}
	    if(symmetric) gsqr.el[j] = fsqr.el[j];
            else dot(G.el[j], G.el[j], &(gsqr.el[j]));
          } 
      }
    
    /* Compute the matrix entries: */
    if (verbose) { fprintf(stderr, "computing matrix...\n"); }
    for (i = 0; i < m; i++)
      { SOFunction *f = F.el[i];
        if (veryverbose) 
          { fprintf(stderr, "[row %d:", i); }
        else if (verbose)
          { fprintf(stderr, "-"); }
        
        for (j = 0; j < n; j++)
          { if (symmetric && (j < i))
              { nat ksym = SOMatrix_Find(Ments.el, kM, j, i);
                if (veryverbose) { fprintf(stderr, " !"); }
                fg = (ksym >= INT_MAX ? 0.0 : Ments.el[ksym].va); 
              }
            else 
              { SOFunction *g = G.el[j];
                if (veryverbose) { fprintf(stderr, " *"); }
                dot(f, g, &fg); //fg = dot(f, g);
              }
            if ((minVal > 0.0) && (fabs(fg) < minVal*sqrt(fsqr.el[i]*gsqr.el[j])))
              { fg = 0.0; }
            if (fg != 0.0)
              { if (veryverbose) { fprintf(stderr, "<F%d,G%d>=%.7f", i, j, fg); }
                MatEntry_vec_expand(&Ments, kM);
                { MatEntry *Mij = &(Ments.el[kM]);
                  Mij->row = i; Mij->col = j;
                  Mij->va = fg; kM++;
                }
              }
          }
        if (veryverbose) { fprintf(stderr, "]\n"); }
      }
    if (veryverbose) { fprintf(stderr, "\n"); }
    MatEntry_vec_trim(&Ments, kM);
    if (minVal > 0.0) { free(fsqr.el); free(gsqr.el); }
    return (SOMatrix){F.nel, G.nel, Ments};
  }
 
void SOBasisMatrices_RecomputeVectorGen
  ( SOFunction* f,
    Basis bas,
    SOFunction_fgDot dot,
    double_vec y,
    bool verbose
  )
  { int i, j, fd = f->d->fDim;
    double sdy2[f->d->fDim], bNew[f->d->fDim], dy;
    for(j = 0; j < fd; j++) sdy2[j] = 0.0;

    if (verbose) { fprintf(stderr, "computing right-hand side \"y\"\n"); }
    for (i = 0; i < bas.nel; i++)
      { dot(f, bas.el[i], bNew);
        for(j = 0; j < fd; j++)
          {
            if (verbose) 
              { dy = bNew[j] - y.el[i*fd + j];
                sdy2[j] += dy*dy;
                fprintf( stderr, "  y[%3d][%3d] = %16.12f  change = %16.12f\n", i, j, bNew[j], dy ); 
              }
            y.el[i*fd + j] = bNew[j];
          }
      }
    if (verbose) 
      for(j = 0; j < fd; j++)
        { fprintf(stderr, "\n");
          fprintf(stderr, "  total change in \"y\"[%3d] = %16.12f\n", j, sqrt(sdy2[j]));
        }
}
   
SOMatrix SOBasisMatrices_BuildMatrixEval
  ( Basis F,
    Basis G, 
    SOFunction *w,
    bool ignoreTriang,
    double minVal,
    bool verbose
  )
  { nat pDim = F.el[0]->d->pDim;
    nat fDim = F.el[0]->d->fDim;
    FuncMap NoMap = IdentityMap(fDim, pDim);
    
    auto void dot(SOFunction *f, SOFunction *g, double *iv);
      /* Computes the dot product of {f} and {g}. */
    
    void dot(SOFunction *f, SOFunction *g, double *iv)
      { if (ignoreTriang)
          { SOFunction_RawDot(f, &NoMap, g, &NoMap, w, NULL, iv); return;}
        else
          { SOFunction_Dot(f, &NoMap, g, &NoMap, w, NULL, iv, FALSE);}
      }
   
    bool symmetric = ((F.nel == G.nel) & (F.el == G.el));
    return SOBasisMatrices_BuildMatrixGen( F, G, dot, minVal, symmetric, verbose);
  }

SOMatrix SOBasisMatrices_BuildMatrixGrad
  ( Basis F, 
    Basis G, 
    SOFunction *w,
    bool ignoreTriang,
    double minVal,
    bool verbose
  )
  { 
    auto void dot(SOFunction *f, SOFunction *g, double *iv);
      /* Computes the dot product of the gradients of {f} and {g}. */
    
    void dot(SOFunction *f, SOFunction *g, double *iv)
      { if (ignoreTriang)
          { SOFunction_RawGradDot(f, g, NULL, iv); return;}
        else
          { SOFunction_GradDot(f, g, w, NULL, iv, FALSE); return;}
      }
    
    bool symmetric = ((F.nel == G.nel) & (F.el == G.el));
    return SOBasisMatrices_BuildMatrixGen( F, G, dot, minVal, symmetric, verbose);
  }

SOMatrix SOBasisMatrices_BuildMatrixGrad1
  ( Basis F, 
    Basis G, 
    SOFunction *w,
    bool ignoreTriang,
    double minVal,
    bool verbose,
    int gradaxis
  )
  { 
    auto void dot(SOFunction *f, SOFunction *g, double *iv);
      /* Computes the dot product of the gradients of {f} and {g}. */
    
    void dot(SOFunction *f, SOFunction *g, double *iv)
      {
/*         if (ignoreTriang) */
/*           { SOFunction_RawGrad1Dot(f, g, NULL, iv, gradaxis); return;} */
/*         else */
          { SOFunction_Grad1Dot(f, g, w, NULL, iv, FALSE, gradaxis); return;}
      }
    
    bool symmetric = FALSE;
    return SOBasisMatrices_BuildMatrixGen( F, G, dot, minVal, symmetric, verbose);
  }
 
void SOBasisMatrices_ChangeBasis
  ( double_vec a,
    SOMatrix FG,
    SOMatrix GGL,
    double_vec b 
  )
  { nat n = FG.cols;
    double_vec c = double_vec_new(n);
    double_vec y = double_vec_new(n);
    SOMatrix_MulRow(a, FG, c);
    SOMatrix_DivCol(GGL, c, y);
    SOMatrix_DivRow(y, GGL, b);
  }