algebra.h

Go to the documentation of this file.

#ifndef ALGEBRAH

#include "basic_algebra.h"

#if !defined(_LAPACK) && !defined(_GSL)
  #error CuikSuite requires LAPACK or GSL
#endif

#define DAMPED_NEWTON 0

#define DN_F 0.1

#define DN_MF 0.95

#define DN_W 0.25

#ifdef _GSL
  #include <gsl/gsl_linalg.h>

  typedef struct {
    unsigned int nr;        
    unsigned int nc;        
    double *Ab;             
    double *bb;             
    gsl_matrix_view A;      
    gsl_vector_view b;      
    gsl_matrix *V;          
    gsl_vector *S;          
    gsl_vector *tmp;        
    gsl_matrix_view VT;     
    gsl_matrix_view UT;     
    gsl_vector_view ST;     
    gsl_vector *w;          
    #if (DAMPED_NEWTON)
      double alpha;         
    #endif
  } TNewton;

  typedef struct {
    unsigned int nr;    
    unsigned int nc;    
    double *Ab;         
    double *bb;         
    double *xb;         
    gsl_matrix_view A;  
    gsl_vector_view b;  
    gsl_vector_view x;  
    gsl_permutation *p; 
    gsl_vector *tau;    
    gsl_vector *res;    
  } TLinearSystem;
#endif

#ifdef _LAPACK
  #ifdef _LAPACKE
    #include <lapacke.h>
  #else
    /* assuming clapack */
    #include <clapack.h>
  #endif

  typedef struct {
    int nr;                 
    int nc;                 
    double *Ab;             
    double *bb;             
    double *s;              
    int lwork;              
    double *work;           
    #if (DAMPED_NEWTON)
      double alpha;         
    #endif
  } TNewton;

  typedef struct {
    int nr;             
    int nc;             
    double *Ab;         
    double *bb;         
    double *xb;         
    int *p;             
    int lwork;          
    double *work;       
  } TLinearSystem;
#endif


/*******************************************************************************/
/* The functions below have different implemenatation depending on the helper  */
/* library being used (gsl, eigen, etc). In other words, these are the         */
/* functions to re-implement to change to another helper library.              */

int MatrixDeterminantSgn(double epsilon,unsigned int n,double *A);

double MatrixDeterminant(unsigned int n,double *A);

unsigned int FindRank(double epsilon,unsigned int nr,unsigned int nc,double *mT);

unsigned int FindKernel(unsigned int nr,unsigned int nc,double *mT,
                        unsigned int dof,boolean check,double epsilon,double **T);

unsigned int FindKernelAndIndependentRows(unsigned int nr,unsigned int nc,double *mT,
                                          unsigned int dof,double epsilon,boolean *singular,
                                          boolean **IR,double **T);

void InitNewton(unsigned int nr,unsigned int nc,TNewton *n);

double *GetNewtonMatrixBuffer(TNewton *n);


double *GetNewtonRHBuffer(TNewton *n);

void NewtonSetMatrix(unsigned int i,unsigned int j,double v,TNewton *n);

void NewtonSetRH(unsigned int i,double v,TNewton *n);


int NewtonStep(double nullSingularValue,double *x,double *dif,
               TNewton *n);

double NewtonGetResult(unsigned int i,TNewton *n);

void DeleteNewton(TNewton *n);

void InitLS(unsigned int nr,unsigned int nc,TLinearSystem *ls);

double *GetLSMatrixBuffer(TLinearSystem *ls);

double *GetLSRHBuffer(TLinearSystem *ls);

double *GetLSSolutionBuffer(TLinearSystem *ls);

void LSSetMatrix(unsigned int i,unsigned int j,double v,TLinearSystem *ls);

void LSSetRH(unsigned int i,double v,TLinearSystem *ls);


int LSSolve(TLinearSystem *ls);

void DeleteLS(TLinearSystem *ls);

#endif

#ifdef _LAPACK



#endif