SampleDrivers/Matrix_Compression_and_Recovery/SMB/sparse_jac_hess.cpp

Go to the documentation of this file.

#include "math.h"
#include "stdlib.h"
#include "stdio.h"

#include "adolc.h"
#include "sparse/sparsedrivers.h"

#define repnum 10

//------------------------------------------------------------------------------------
// for time measurements

#include <sys/time.h>
#include "ColPackHeaders.h"

using namespace ColPack;

double k_getTime() {
   struct timeval v;
   struct timezone z;
   gettimeofday(&v, &z);
   return ((double)v.tv_sec)+((double)v.tv_usec)/1000000;
}

//------------------------------------------------------------------------------------
// required for second method

using namespace std;

#include <list>
#include <map>
#include <string>
#include <vector>

//------------------------------------------------------------------------------------
// as before

#define tag_f 1
#define tag_red 2
#define tag_HP 3

#define tag_c 4

// problem definition -> eval_fun.c
void init_dim(int *n, int *m);
void init_startpoint(double *x, int n);
double  feval(double *x, int n);
adouble feval_ad(double *x, int n);
void ceval(double *x, double *c, int n);
void ceval_ad(double *x, adouble *c, int n);
adouble feval_ad_mod(double *x, int n);
adouble feval_ad_modHP(double *x, int n);

void printmat(char* kette, int n, int m, double** M);
void printmatint(char* kette, int n, int m, int** M);
void printmatint_c(char* kette, int m,unsigned int** M);


int main()
{
  int i, j, k, l, sum, n, m, nnz, direct = 1, found;
  double f;
  double *x, *c;
  adouble fad, *xad, *cad;
  //double** Zpp;
  //double** Zppf;
  double** J;
  //double* s;
  //int p_H_dir, p_H_indir;
  int tape_stats[11];
  int num;
  FILE *fp_JP;

  double **Seed_J;
  double **Jc;
  int p_J;

  int recover = 1;
  int jac_vec = 1;
  int compute_full = 1;
  int output_sparsity_pattern_J = 0;
  //int output_sparsity_pattern_H = 1;
  //int use_direct_method = 1;
  //int output_direct = 0;
  //int use_indirect_method = 1;
  //int output_indirect = 0;

  double t_f_1, t_f_2, div_f=0, div_c=0, div_JP=0, div_JP2=0, div_Seed=0, div_Seed_C=0, div_Jc=0, div_Jc_C=0, div_rec=0, div_rec_C=0, div_J=0;
  //double test;
  unsigned int *rind;
  unsigned int *cind;
  double *values;

  //tring s_InputFile = "test_mat.mtx";
  //string s_InputFile = "jac_pat.mtx";


//------------------------------------------------------------------------------------
// problem definition + evaluation

  init_dim(&n,&m); // initialize n and m

  printf(" n = %d m = %d\n",n,m);

  x =   (double*)  malloc(n*sizeof(double)); // x: vector input for function evaluation
  c =   (double*)  malloc(m*sizeof(double)); // c: constraint vector
  cad = new adouble[m];

  init_startpoint(x,n);

  t_f_1 = k_getTime();
  for(i=0;i<repnum;i++)
    f = feval(x,n);
  t_f_2 = k_getTime();
  div_f = (t_f_2 - t_f_1)*1.0/repnum;
  printf("XXX The time needed for function evaluation:  %10.6f \n \n", div_f);


  t_f_1 = k_getTime();
  for(i=0;i<repnum;i++)
    ceval(x,c,n);
  t_f_2 = k_getTime();
  div_c = (t_f_2 - t_f_1)*1.0/repnum;
  printf("XXX The time needed for constraint evaluation:  %10.6f \n \n", div_c);


  trace_on(tag_f);

    fad = feval_ad(x, n); // feval_ad: derivative of feval

    fad >>= f;

  trace_off();

  trace_on(tag_c);

    ceval_ad(x, cad, n); //ceval_ad: derivative of ceval

    for(i=0;i<m;i++)
      cad[i] >>= f;

  trace_off();
  //return 1;

  tapestats(tag_c,tape_stats);              // reading of tape statistics
  printf("\n    independents   %d\n",tape_stats[0]);
  printf("    dependents     %d\n",tape_stats[1]);
  printf("    operations     %d\n",tape_stats[5]);
  printf("    buffer size    %d\n",tape_stats[4]);
  printf("    maxlive        %d\n",tape_stats[2]);
  printf("    valstack size  %d\n\n",tape_stats[3]);


//------------------------------------------------------------------------------------
// full Jacobian:

  div_J = -1;

  if(compute_full == 1)
    {
      J =  myalloc2(m,n);

      t_f_1 = k_getTime();
      jacobian(tag_c,m,n,x,J);
      t_f_2 = k_getTime();
      div_J = (t_f_2 - t_f_1);

      printf("XXX The time needed for full Jacobian:  %10.6f \n \n", div_J);
      printf("XXX runtime ratio:  %10.6f \n \n", div_J/div_c);

      //save the matrix into a file (non-zero entries only)

        fp_JP = fopen("jac_full.mtx","w");

        fprintf(fp_JP,"%d %d\n",m,n);

        for (i=0;i<m;i++)
          {
            for (j=0;j<n;j++)
              if(J[i][j]!=0.0) fprintf(fp_JP,"%d %d %10.6f\n",i,j,J[i][j] );
          }
        fclose(fp_JP);
    }

//------------------------------------------------------------------------------------
  printf("XXX THE 4 STEP TO COMPUTE SPARSE MATRICES USING ColPack \n \n");
// STEP 1: Determination of sparsity pattern of Jacobian JP:

  unsigned int  *rb=NULL;          /* dependent variables          */
  unsigned int  *cb=NULL;          /* independent variables        */
  unsigned int  **JP=NULL;         /* compressed block row storage */
  int ctrl[2];

  JP = (unsigned int **) malloc(m*sizeof(unsigned int*));
  ctrl[0] = 0; ctrl[1] = 0;


  t_f_1 = k_getTime();
  jac_pat(tag_c, m, n, x, JP, ctrl);    //ADOL-C calculate the sparsity pattern
  t_f_2 = k_getTime();
  div_JP = (t_f_2 - t_f_1);

  printf("XXX STEP 1: The time needed for Jacobian pattern:  %10.6f \n \n", div_JP);
  printf("XXX STEP 1: runtime ratio:  %10.6f \n \n", div_JP/div_c);


  nnz = 0;
  for (i=0;i<m;i++)
    nnz += JP[i][0];

  printf(" nnz %d \n",nnz);
  printf(" hier 1a\n");


//------------------------------------------------------------------------------------
// STEP 2: Determination of Seed matrix:

  double tg_C;
  int dummy;

  t_f_1 = k_getTime();

  BipartiteGraphPartialColoringInterface * gGraph = new BipartiteGraphPartialColoringInterface(SRC_MEM_ADOLC, JP, m, n);
  //gGraph->PrintBipartiteGraph();
  t_f_2 = k_getTime();

  printf("XXX STEP 2: The time needed for Graph construction:  %10.6f \n \n", (t_f_2-t_f_1) );
  printf("XXX STEP 2: runtime ratio:  %10.6f \n \n", (t_f_2-t_f_1)/div_c);
  
  t_f_1 = k_getTime();
  //gGraph->GenerateSeedJacobian(&Seed_J, &dummy, &p_J,
  //                          "NATURAL", "COLUMN_PARTIAL_DISTANCE_TWO");
  gGraph->PartialDistanceTwoColoring("NATURAL", "COLUMN_PARTIAL_DISTANCE_TWO");
  t_f_2 = k_getTime();

  printf("XXX STEP 2: The time needed for Coloring:  %10.6f \n \n", (t_f_2-t_f_1));
  printf("XXX STEP 2: runtime ratio:  %10.6f \n \n", (t_f_2-t_f_1)/div_c);
  
  t_f_1 = k_getTime();
  Seed_J = gGraph->GetSeedMatrix(&dummy, &p_J);
  t_f_2 = k_getTime();
  tg_C = t_f_2 - t_f_1;


  printf("XXX STEP 2: The time needed for Seed generation:  %10.6f \n \n", tg_C);
  printf("XXX STEP 2: runtime ratio:  %10.6f \n \n", tg_C/div_c);

  //*/

//------------------------------------------------------------------------------------
// STEP 3: Jacobian-matrix product:

// ADOL-C:
//*
  if (jac_vec == 1)
    {

      Jc = myalloc2(m,p_J);
      t_f_1 = k_getTime();
      printf(" hier 1\n");
      fov_forward(tag_c,m,n,p_J,x,Seed_J,c,Jc);
      printf(" hier 2\n");
      t_f_2 = k_getTime();
      div_Jc = (t_f_2 - t_f_1);

      printf("XXX STEP 3: The time needed for Jacobian-matrix product:  %10.6f \n \n", div_Jc);
      printf("XXX STEP 3: runtime ratio:  %10.6f \n \n", div_Jc/div_c);


    }

//------------------------------------------------------------------------------------
// STEP 4: computed Jacobians/ recovery


  if (recover == 1)
    {


      JacobianRecovery1D jr1d;

      printf("m = %d, n = %d, p_J = %d \n",m,n,p_J);
      //printmatint_c("JP Jacobian Pattern",m,JP);
      //printmat("Jc Jacobian compressed",m,p_J,Jc);

      t_f_1 = k_getTime();
      jr1d.RecoverD2Cln_CoordinateFormat (gGraph, Jc, JP, &rind, &cind, &values);
      t_f_2 = k_getTime();
      div_rec_C = (t_f_2 - t_f_1);

      printf("XXX STEP 4: The time needed for Recovery:  %10.6f \n \n", div_rec_C);
      printf("XXX STEP 4: runtime ratio:  %10.6f \n \n", div_rec_C/div_c);

      //save recovered matrix into file

        fp_JP = fopen("jac_recovered.mtx","w");

        fprintf(fp_JP,"%d %d %d\n",m,n,nnz);

        for (i=0;i<nnz;i++)
          {
              fprintf(fp_JP,"%d %d %10.6f\n",rind[i],cind[i],values[i] );
          }
        fclose(fp_JP);

    }

    /*By this time, if you compare the 2 output files: jac_full.mtx and jac_recovered.mtx
    You should be able to see that the non-zero entries are identical
    */


  free(JP);
  delete[] cad;
  free(c);
  free(x);
  delete gGraph;
  if(jac_vec == 1) {
    myfree2(Jc);
  }
  if(compute_full == 1)
    {
      myfree2(J);
    }

  return 0;

}




/***************************************************************************/

void printmat(char* kette, int n, int m, double** M)
{ int i,j;

  printf("%s \n",kette);
  for(i=0; i<n ;i++)
  {
    printf("\n %d: ",i);
    for(j=0;j<m ;j++)
        printf(" %10.4f ", M[i][j]);
  }
  printf("\n");
}

void printmatint(char* kette, int n, int m, int** M)
{ int i,j;

  printf("%s \n",kette);
  for(i=0; i<n ;i++)
  {
    printf("\n %d: ",i);
    for(j=0;j<m ;j++)
      printf(" %d ", M[i][j]);
  }
  printf("\n");
}


void printmatint_c(char* kette, int m,unsigned int** M)
{ int i,j;

  printf("%s \n",kette);
  for (i=0;i<m;i++)
    {
    printf("\n %d (%d): ",i,M[i][0]);
      for (j=1;j<=M[i][0];j++)
        printf("\t%d ",M[i][j]);
    }
  printf("\n");
}