#include "tao.h"
static  char help[]="t-logistic regression using TAO \n\
Input parameters are:\n\
  -data    <filename> : file which contains data in Petsc binary format\n\\n";

#undef __FUNCT__
#define __FUNCT__ "main"
int main(int argc,char **argv)
{
  PetscErrorCode     info;              
  char               data_path[PETSC_MAX_PATH_LEN];
  PetscTruth         data_flg;
  PetscViewer        view;
  Vec                w;                     // solution vector
  Mat data,data_parts;
  Vec fx;
  PetscInt m,dim,nclasses;
  PetscInt m_local,dim_local;
  
  nclasses=3;
  
  // Initialize TAO and PETSc
  PetscInitialize(&argc,&argv,(char *)0,help);
  
  // Check for command line arguments to override defaults
  info=PetscOptionsGetString(PETSC_NULL,"-data",data_path,PETSC_MAX_PATH_LEN-1,&data_flg);CHKERRQ(info);
  
  PetscViewerBinaryOpen(PETSC_COMM_WORLD,data_path,FILE_MODE_READ,&view);
  MatLoad(view,MATAIJ,&data_parts);
  PetscViewerDestroy(view);
  
  info=MatCreateMAIJ(data_parts,nclasses,&data);CHKERRQ(info);
  info=MatGetSize(data_parts,&m,&dim);CHKERRQ(info);
  info=MatGetLocalSize(data,&m_local,&dim_local);CHKERRQ(info);

  // Allocate space for various objects
  info=VecCreate(PETSC_COMM_WORLD,&w);CHKERRQ(info);
  info=VecSetSizes(w,dim_local,dim*nclasses);CHKERRQ(info);
  info=VecSetFromOptions(w);CHKERRQ(info);
  info=PetscObjectSetName((PetscObject) w,"Solution");CHKERRQ(info);
  
  for (PetscInt i=0; i<nclasses*dim; i++) {
    PetscScalar v=i;
    info = VecSetValues(w,1,&i,&v,INSERT_VALUES);CHKERRQ(info);
  }

  info=VecAssemblyBegin(w);CHKERRQ(info);
  info=VecAssemblyEnd(w);CHKERRQ(info);
  
  info=VecCreate(PETSC_COMM_WORLD,&fx);CHKERRQ(info);
  info=VecSetSizes(fx,m_local,m*nclasses);CHKERRQ(info);
  info=VecSetFromOptions(fx);CHKERRQ(info);
  info=PetscObjectSetName((PetscObject) fx,"function");CHKERRQ(info);


  // info=VecView(w,PETSC_VIEWER_STDOUT_WORLD);CHKERRQ(info);
  info=MatGetSize(data,&m,&dim);CHKERRQ(info);
  PetscPrintf(PETSC_COMM_WORLD,"Number of data points: %d \n",m);  
  PetscPrintf(PETSC_COMM_WORLD,"Dimensions: %d \n",dim);
  
  info=MatGetLocalSize(data,&m_local,&dim_local);CHKERRQ(info);
  PetscSynchronizedPrintf(PETSC_COMM_WORLD,"local m: %d local dim: %d \n",m_local, dim_local);  
  PetscSynchronizedFlush(PETSC_COMM_WORLD);


  info=MatMult(data,w,fx);CHKERRQ(info);
  // info=VecView(fx,PETSC_VIEWER_STDOUT_WORLD);CHKERRQ(info);
  
  // Free PETSc data structures 
  info=VecDestroy(w);CHKERRQ(info);
  info=VecDestroy(fx);CHKERRQ(info);
  info=MatDestroy(data_parts);CHKERRQ(info);
  info=MatDestroy(data);CHKERRQ(info);

  // Finalize Petsc
  PetscFinalize();

  return 0;
}