jacobian.c

Go to the documentation of this file.

#include "jacobian.h"

#include "error.h"
#include "basic_algebra.h"

void InitJacobian(Tvariables *vs,Tequations *eqs,TJacobian *j)
{
  unsigned int i;

  j->neqs=NEqualityEquations(eqs);
  j->nvars=NVariables(vs);

  if (j->nvars==0)
    Error("Empty problem in InitJacobian");

  NEW(j->J,j->nvars,Tequations);
  for(i=0;i<j->nvars;i++)
    {
      DeriveEqualityEquations(i,&(j->J[i]),eqs);
      CacheScalarEQUInfo(&(j->J[i]));
    }
}

void CopyJacobian(TJacobian *j_dst,TJacobian *j_src)
{
  unsigned int i;

  j_dst->neqs=j_src->neqs;
  j_dst->nvars=j_src->nvars;

  NEW(j_dst->J,j_dst->nvars,Tequations);
  for(i=0;i<j_dst->nvars;i++)
    {
      CopyEquations(&(j_dst->J[i]),&(j_src->J[i]));
      CacheScalarEQUInfo(&(j_dst->J[i]));
    }
}

void GetJacobianSize(unsigned int *nr,unsigned int *nc,TJacobian *j)
{
  *nr=j->neqs;
  *nc=j->nvars;
}

Tequations *GetJacobianColumn(unsigned int nv,TJacobian *j)
{
  if (nv>=j->nvars)
    Error("Index out of range in GetJacobianColumn");
  
  return(&(j->J[nv]));
}

Tequation *GetJacobianEquation(unsigned int r,unsigned int c,TJacobian *j)
{
  if ((r>=j->neqs)||(c>=j->nvars))
    Error("Index out of range in GetJacobianEquation");

  return(GetEquation(r,&(j->J[c])));
}

void AllocateJacobianEvaluation(double ***m,TJacobian *j)
{
  if (j->neqs>0)
    {
      unsigned int i;
  
      NEW(*m,j->neqs,double *);
      for(i=0;i<j->neqs;i++)
        {
          NEW((*m)[i],j->nvars,double);
        }
    }
}

void EvaluateJacobian(double *v,double **m,TJacobian *j)
{
  if (j->neqs>0)
    {
      unsigned int i,k;
      double *e;

      NEW(e,j->neqs,double);
      for(i=0;i<j->nvars;i++)
        {
          EvaluateEqualityEquations(FALSE,v,e,&(j->J[i]));
          for(k=0;k<j->neqs;k++)
            m[k][i]=e[k];
        }
      free(e);
    }
}

void EvaluateJacobianInVector(double *v,unsigned int nr,unsigned int nc,double *m,TJacobian *j)
{
  if (j->neqs>0)
    {
      unsigned int i;
      double *e;

      NEW(e,j->neqs,double);
      for(i=0;i<j->nvars;i++)
        {
          EvaluateEqualitySparseEquations(v,e,&(j->J[i]));
          SubMatrixFromMatrix(j->neqs,1,e,0,i,nr,nc,m);
        }
      free(e);
    }
}

void EvaluateJacobianSubSetInVector(double *v,boolean *sr,unsigned int nr,unsigned int nc,double *m,TJacobian *j)
{
  if (j->neqs>0)
    {
      unsigned int i,l;
      double *e;
  
      l=0;
      for(i=0;i<j->neqs;i++)
        {
          if (sr[i]) l++;
        }

      NEW(e,l,double);
      for(i=0;i<j->nvars;i++)
        {
          EvaluateSubSetEqualitySparseEquations(v,sr,e,&(j->J[i]));
          SubMatrixFromMatrix(l,1,e,0,i,nr,nc,m);
        }
      free(e);
    }
}

void EvaluateTransposedJacobianInVector(double *v,unsigned int nr,unsigned int nc,double *m,TJacobian *j)
{
  if (j->neqs>0)
    {
      unsigned int i;
      double *e;

      NEW(e,j->neqs,double);
      for(i=0;i<j->nvars;i++)
        {
          EvaluateEqualitySparseEquations(v,e,&(j->J[i]));
          SubMatrixFromMatrix(1,j->neqs,e,i,0,nr,nc,m);
        }
      free(e);
    }
}

void EvaluateTransposedJacobianSubSetInVector(double *v,boolean *sr,unsigned int nr,unsigned int nc,double *m,TJacobian *j)
{
  if (j->neqs>0)
    {
      unsigned int i,l;
      double *e;

      l=0;
      for(i=0;i<j->neqs;i++)
        {
          if (sr[i]) l++;
        }

      NEW(e,l,double);
      for(i=0;i<j->nvars;i++)
        {
          EvaluateSubSetEqualitySparseEquations(v,sr,e,&(j->J[i]));
          SubMatrixFromMatrix(1,l,e,i,0,nr,nc,m);
        }
      free(e);
    }
}

void PrintJacobianEvaluation(FILE *f,double **m,TJacobian *j)
{
  unsigned int i,k;

  fprintf(f,"[");
  for(i=0;i<j->neqs;i++)
    {
      for(k=0;k<j->nvars;k++)
        {
          fprintf(f,"%.16f ",m[i][k]);
        }
      fprintf(f,"\n");
    }
  fprintf(f,"];\n");
}


void FreeJacobianEvaluation(double **m,TJacobian *j)
{
  if (j->neqs>0)
    {
      unsigned int i;

      for(i=0;i<j->neqs;i++)
        free(m[i]);
      free(m);
    }
}

void EvaluateJacobianXVectors(double *p,
                              unsigned int n,unsigned int ng,unsigned int *g,double *v,
                              unsigned int *nr,unsigned int *nc,double ***m,
                              TJacobian *j)
{
  unsigned int i,k;
  double *e;

  /* Size of the output matrix */
  *nr=3*n;
  *nc=j->nvars;

  /* Allocate the output matrix. */
  NEW(*m,3*n,double*);
  for(i=0;i<*nr;i++)
    { NEW((*m)[i],*nc,double); }

  /* Space for the evaluation of each column */
  NEW(e,*nr,double);

  /* Evaluate column-wise */
  for(i=0;i<j->nvars;i++)
    {
      EvaluateEquationsXVectors(p,ng,g,v,e,&(j->J[i]));
      for(k=0;k<(*nr);k++)
        (*m)[k][i]=e[k];
    }

  /* Release space */
  free(e);
}

void PrintJacobian(FILE *f,char **varNames,TJacobian *j)
{
  unsigned int i;

  for(i=0;i<j->nvars;i++)
    {
      fprintf(f,"Jacobian column %u\n",i);
      PrintEquations(f,varNames,&(j->J[i]));
    }
}

void DeleteJacobian(TJacobian *j)
{
  unsigned int i;

  for(i=0;i<j->nvars;i++)
    DeleteEquations(&(j->J[i]));
  free(j->J);
}