MatrixBlockIP.h

Go to the documentation of this file.

#ifndef MATRIXBLOCKIP_H
#define MATRIXBLOCKIP_H

#include "SparseChol.h"
#include <vector>
#include <fstream>

using namespace std;

class MatrixBlockIP {
 public:
  
  /* 
     For every type of matrix we need routines to:
     - multiply matrix by vector, and matrix transposed by vector
     - solve systems with matrix (either through Cholesky or specialized routines)
     - routine to transform matrix to a GENERAL matrix (to get full matrix of constraints)
  */

  TYPE_MATRIX type; 
  int m; 
  int n; 
  bool ija; 
  bool rowwise; 
  bool columnwise; 
  
  int nz; 
  double *a; 
  
  int *inirowa; 
  int *icola; 
  
  int *inicola; 
  int *irowa; 
  
  TYPE_ORIENTATION type_orientation; 
  int num_arcs; 
  int num_nodes; 
  int *src; 
  int *dst; 
  
  double *d1; 
  double *d2; 

  SparseChol chol; 
  
  int sizeL;
  int *inirowLLt;
  // The following vectors are allocated with ALLOC and freed with FREE
  int *icolLLt;
  int *blockD;
  int *inivalD;
  int *icolD;
  double *valD;

  bool made_symbfct_MMt; // for symbolic factorization of M*Theta*M'
  bool made_symbfct_M;   // for symbolic factorization of M
  bool made_analyze_D; 

  int num_blocks; 

  // Constructor
  MatrixBlockIP(int blocks=1);
  // Copy constructor
  MatrixBlockIP(MatrixBlockIP *mbip);
  // Destructor
  ~MatrixBlockIP();
  
  // Copy
  void copy(MatrixBlockIP *mbip);
  
  // Comes back to the initial state
  void reset();
  
  // Comes back to the state after create the matrix
  void restore();
  
  /* Native methods to create different matrix types */
  // Native method to create a general row-wise packed matrix
  void create_general_matrix_row_wise(int m, int n, int nz, int* &inirowa, int* &icola, double* &a);
  // Native method to create a general column-wise packed matrix
  void create_general_matrix_column_wise(int m, int n, int nz, int* &inicola, int* &irowa, double* &a);
  // Native method to create a network matrix
  void create_network_matrix(int num_arcs, int num_nodes, int* &src, int* &dst, bool oriented= true);
  // Native method to create a identity matrix
  void create_identity_matrix(int dim);
  // Native method to create a identity-identity matrix
  void create_idty_idty_matrix(int nrows);
  // Native method to create a diagonal matrix
  void create_diagonal_matrix(int dim, double* &d);
  // Native method to create a diagonal-diagonal matrix
  void create_diag_diag_matrix(int nrows, double* &d1, double* &d2);
/*  
  // Native method to load a general row-wise packed matrix
  void load_general_matrix_row_wise(int m, int n, int nz, int inirowa[], int icola[], double a[], bool ordered=false, bool copy_vectors=true);
  // Native method to load a general column-wise packed matrix
  void load_general_matrix_column_wise(int m, int n, int nz, int inicola[], int irowa[], double a[], bool ordered=false, bool copy_vectors=true);
  // Native method to load a network matrix
  void load_network_matrix(int num_arcs, int num_nodes, int src[], int dst[], bool oriented= true, bool copy_vectors=true);
  // Native method to load a identity matrix
  void load_identity_matrix(int dim);
  // Native method to load a identity-identity matrix
  void load_idty_idty_matrix(int dim);
  // Native method to load a diagonal matrix
  void load_diagonal_matrix(int dim, double d[], bool copy_vectors=true);
  // Native method to load a diagonal-diagonal matrix
  void load_diag_diag_matrix(int dim, double d1[], double d2[], bool copy_vectors=true);
*/
    
  /* Alternative methods to load/build matrices */
  // Creates a general matrix in format ija
  void create_general_matrix_format_ija(int m, int n, int nz, int* &row_index, int* &col_index, double* &values);
/*
  // Loads a general matrix in format ija
  void load_general_matrix_format_ija(int m, int n, int nz, int row_index[]= NULL, int col_index[]=NULL, double values[]= NULL, bool copy_vectors=true);  
*/
  // Builds a packed rowwise format structured matrix based on diagonal blocks and linking constraints blocks
  void compute_full_matrix(int numBlocks, bool sameN, MatrixBlockIP N[], bool sameL, MatrixBlockIP L[], int numActiveLnk=0, int *listActiveLnk=NULL);
  // Create internal structure arrays to compute D = Theta_(numBlocks+1) + sum{i..numBlocks} L_i*Theta_i*L_i'
  void analyze_D(int numBlocks, bool sameL, MatrixBlockIP L[]);
  void analyze_D_diagonal(int numBlocks, bool sameL, MatrixBlockIP L[]);
  void analyze_D_gen_sym_uptr(int numBlocks, bool sameL, MatrixBlockIP L[]);
  // Compute D = Theta_(numBlocks+1) + sum{i..numBlocks} L_i*Theta_i*L_i'
  void compute_D(int numBlocks, bool sameL, MatrixBlockIP L[], bool isActive[], int iniTheta[], double Theta[]);
  void compute_D_diagonal(int numBlocks, bool sameL, MatrixBlockIP L[], bool isActive[], int iniTheta[], double Theta[]);
  void compute_D_gen_sym_uptr(int numBlocks, bool sameL, MatrixBlockIP L[], bool isActive[], int iniTheta[], double Theta[]);
  
    
  /* Matrix transformations */
  // Calculates and stores the matrix in packed rowwise format
  void native_to_general(bool delete_native_format=false);
  // Calculates and stores the matrix in packed rowwise format
  void ija_to_rowwise();
  // Calculates and stores the network matrix (either oriented or nonoriented) in packed rowwise format
  void network_to_general();
  // Calculates and stores the I matrix in packed rowwise format
  void identity_to_general();
  // Calculates and stores the diagonal D matrix in packed rowwise format
  void diagonal_to_general();
  // Calculates and stores the matrix [I I] in packed rowwise format
  void idty_idty_to_general();
  // Calculates and stores the matrix [D1 D2] in packed rowwise format
  void diag_diag_to_general();
  
  // Calculates and stores the network matrix in ija format.
  void network_to_ija_format();
  // Calculates and stores the identity I matrix in ija format.
  void identity_to_ija_format();
  // Calculates and stores the diagonal D1 matrix in ija format.
  void diagonal_to_ija_format();
  // Calculates and stores the matrix  [I I] in ija format
  void idty_idty_to_ija_format();
  // Calculates and stores the matrix  [D1 D2] in ija format
  void diag_diag_to_ija_format();
  // Order the matrix when has been created in row-wise or column-wise format
  void order_matrix();
  
  /* Matrix operations */
  // Matrix-vector product vout=M*vin (vout(m),vin(n))
  void mul_Mv(double vout[], const double vin[]);
  // MatrixTranspose-vector product vout=M(t)*vin (vout(n),vin(m)
  void mul_Mtv(double vout[], const double vin[]);
  // Add matrix-vector product vout += M*vin (vout(m),vin(n))
  void add_mul_Mv(double vout[], const double vin[]);
  // Add matrixTranspose-vector product vout += M(t)*vin (vout(n),vin(m)
  void add_mul_Mtv(double vout[], const double vin[]);
  // Given a column determines if exists an element for each row of the matrix, in case of exist returns the value.
  void exist_var_in_row(int column, bool appear[], double coefs[]);
  // Adds a new column into the matrix
  void add_new_column(int size, int irows[], double values[]);
  // Adds new columns into the matrix
  void add_new_columns(int num_columns, int size[], int *irows[], double *values[]);
  // Changes the sign of some columns
  void change_columns_sign(int size, int columns[]);
  // Changes the sign of some rows
  void change_rows_sign(int size, int rows[]);
  // Deletes some rows
  void delete_rows(int size, int rows[]);
  // Gets a column of the matrix
  int get_column(int column, int* &irows, double* &values, bool invert_sign= false);
  
  // SparseChol interface functions
  // next three for systems with M*Theta*M'
  void symbolic_fact_MMt(CHOL_SOLVER chslv=SPRSBLKLLT, int *prov_pfa= NULL, int *prov_ipfa= NULL, NUMBERING prov_pfa_numbering= NOT_COMPUTED);
  void numeric_fact_MMt(double *Theta, int i_k=0);
  void numeric_solve_MMt(double *rhs, int i_k= 0, WHO_PERMUTES whoperm= CHOLESKY);
  // next three for systems with M (M positive semidefinite matrix)
  void symbolic_fact_M(CHOL_SOLVER chslv=SPRSBLKLLT, int *prov_pfa= NULL, int *prov_ipfa= NULL, NUMBERING prov_pfa_numbering= NOT_COMPUTED);
  void numeric_fact_M();
  void numeric_solve_M(double *rhs, WHO_PERMUTES whoperm= CHOLESKY);

  // access to data in class
  int get_pfa(int i); // returns SparseChol pfa[i] with no check
  int get_ipfa(int i); // returns SparseChol ipfa[i] with no check
  int get_maxlnz(); // returns SparseChol maxlnz with no check
  int get_maxfillin(); // returns SparseChol maxfillin with no check
  int get_njka(); // returns SparseChol njka with no check
  int get_num_zero_pivots(); // returns SparseChol num_zero_pivots with no check
  int get_num_semidef_matrix(); // returns SparseChol num_semidef_matrix with no check
  
  // only to test
  void print_matrix(bool ija= true, bool start_one= true);
  void print_matrix(ofstream &outfile, bool print_ija= true, bool start_one= true);
  // Print some positions and name of a vector
  void print_vector(double *v, int size, string name);

  // Pre: icola is ordered
  // Convert a matrix in column-wise format to row-wise format.
  void column_wise_to_row_wise_format();
  
  // Pre: irowa is ordered
  // Convert a matrix in row-wise format to column-wise format.
  void row_wise_to_column_wise_format();
    
  private:
  
  // free all the memory allocated by the application - not the user
  void free_memory();

  void mul_Mv_row_wise(double vout[], const double vin[]);
  void mul_Mtv_row_wise(double vout[], const double vin[]);
  void add_mul_Mv_row_wise(double vout[], const double vin[]);
  void add_mul_Mtv_row_wise(double vout[], const double vin[]);
  
  void mul_Mv_column_wise(double vout[], const double vin[]);
  void mul_Mtv_column_wise(double vout[], const double vin[]);
  void add_mul_Mv_column_wise(double vout[], const double vin[]);
  void add_mul_Mtv_column_wise(double vout[], const double vin[]);

  void mul_Mv_network(double vout[], const double vin[]);
  void mul_Mtv_network(double vout[], const double vin[]);
  void add_mul_Mv_network(double vout[], const double vin[]);
  void add_mul_Mtv_network(double vout[], const double vin[]);

  void mul_Mv_identity(double vout[], const double vin[]);
  void add_mul_Mv_identity(double vout[], const double vin[]);

  void mul_Mv_diagonal(double vout[], const double vin[]);
  void add_mul_Mv_diagonal(double vout[], const double vin[]);

  void mul_Mv_idty_idty(double vout[], const double vin[]);
  void mul_Mtv_idty_idty(double vout[], const double vin[]);
  void add_mul_Mv_idty_idty(double vout[], const double vin[]);
  void add_mul_Mtv_idty_idty(double vout[], const double vin[]);

  void mul_Mv_diag_diag(double vout[], const double vin[]);
  void mul_Mtv_diag_diag(double vout[], const double vin[]);
  void add_mul_Mv_diag_diag(double vout[], const double vin[]);
  void add_mul_Mtv_diag_diag(double vout[], const double vin[]);
  
  void mul_Mv_gen_sym_uptr(double vout[], const double vin[]); // M=M', then mul_Mtv call this one
  void add_mul_Mv_gen_sym_uptr(double vout[], const double vin[]); // M=M', then mul_Mtv call this one

  // Order a matrix packed format (both column-wise and row-wise)
  void order_packed_format(int begsize, int nz, int* &beg, int* &ind, double* &val);
  
  struct Order_ija {
    int *row; 
    int *col; 
    bool operator() (int i,int j) 
    {
      return ((row[i] < row[j]) || ((row[i] == row[j]) && (col[i] < col[j])));
    }
  };
  
  struct Order_vector {
    int *v; 
    bool operator() (int i,int j) 
    {
      return (v[i] < v[j]);
    }
  };
};


inline void MatrixBlockIP::symbolic_fact_MMt(CHOL_SOLVER chslv, int *prov_pfa, int *prov_ipfa, NUMBERING prov_pfa_numbering)
{
  chol.initialize(chslv,type,type_orientation);

  if (m <= 0) return; // for 0 rows matrices do nothing

  switch (type) {
  case GENERAL:
    chol.symbolic_fact_MMt(m, n, nz, icola, inirowa, a, num_blocks, prov_pfa, prov_ipfa, prov_pfa_numbering);
    break;
  case NETWORK:
    chol.symbolic_fact_MMt(num_nodes, num_arcs, src, dst, num_blocks, prov_pfa, prov_ipfa, prov_pfa_numbering);
    break;
  case IDENTITY:
  case IDTY_IDTY:
    chol.symbolic_fact_MMt(m, num_blocks);
    break;
  case DIAGONAL:
    chol.symbolic_fact_MMt(m,d1,num_blocks);
    break;
  case DIAG_DIAG:
    chol.symbolic_fact_MMt(m,d1,d2,num_blocks);
    break;
  default:
    throw ExceptionBlockIP(NOT_FOR_THIS_MATRIX_TYPE, __FILE__, __LINE__);
  }
}

inline void MatrixBlockIP::numeric_fact_MMt(double *Theta, int i_k)
{
  if (m <= 0) return; // for 0 rows matrices do nothing

  // if first time then creation and symbolic factorization of M*Theta*M'
  if (!made_symbfct_MMt) {
    symbolic_fact_MMt();
    made_symbfct_MMt= true;
  }

  chol.numeric_fact_MMt(Theta, i_k);
}

inline void MatrixBlockIP::numeric_solve_MMt(double *rhs, int i_k, WHO_PERMUTES whoperm)
{
  if (m <= 0) return; // for 0 rows matrices do nothing

  chol.numeric_solve_MMt(rhs, i_k, whoperm);
}

inline void MatrixBlockIP::symbolic_fact_M(CHOL_SOLVER chslv, int *prov_pfa, int *prov_ipfa, NUMBERING prov_pfa_numbering)
{
  chol.initialize(chslv,type,type_orientation);

  if (m <= 0) return; // for 0 rows matrices do nothing

  // only implemented (and needed by now) for GEN_SYM_UPTR and DIAGONAL matrices
  switch (type) {
  case DIAGONAL:
    chol.symbolic_fact_M(m);
    break;
  case GEN_SYM_UPTR:
    chol.symbolic_fact_M(m, nz, icola, inirowa, prov_pfa, prov_ipfa, prov_pfa_numbering);
    break;
  default:
    throw ExceptionBlockIP(NOT_FOR_THIS_MATRIX_TYPE, __FILE__, __LINE__);
  }
}

inline void MatrixBlockIP::numeric_fact_M()
{
  if (m <= 0) return; // for 0 rows matrices do nothing

  // if first time then symbolic factorization of M
  if (!made_symbfct_M) {
    symbolic_fact_M();
    made_symbfct_M= true;
  }

  switch (type) {
  case DIAGONAL:
    chol.numeric_fact_M(d1);
    break;
  case GEN_SYM_UPTR:
    chol.numeric_fact_M(a);
    break;
  default:
    throw ExceptionBlockIP(NOT_FOR_THIS_MATRIX_TYPE, __FILE__, __LINE__);
  }
}

inline void MatrixBlockIP::numeric_solve_M(double *rhs, WHO_PERMUTES whoperm)
{
  if (m <= 0) return; // for 0 rows matrices do nothing

  chol.numeric_solve_M(rhs, whoperm);
}

inline int MatrixBlockIP::get_pfa(int i) 


{
  return(chol.get_pfa(i));
}

inline int MatrixBlockIP::get_ipfa(int i)

{
  return(chol.get_ipfa(i));
}

inline int MatrixBlockIP::get_maxlnz()

{
  return (chol.get_maxlnz());
}

inline int MatrixBlockIP::get_maxfillin()

{
  return (chol.get_maxfillin());
}

inline int MatrixBlockIP::get_njka()

{
  return (chol.get_njka());
}

inline int MatrixBlockIP::get_num_zero_pivots()

{
  return (chol.get_num_zero_pivots());
}

inline int MatrixBlockIP::get_num_semidef_matrix()

{
  return (chol.get_num_semidef_matrix());
}

inline void MatrixBlockIP::native_to_general(bool delete_native_format)

{
  if (rowwise) return; // already stored

  switch (type) {
  case GENERAL:
    if (columnwise)
      column_wise_to_row_wise_format();
    else if (ija)
      ija_to_rowwise();
    else {
      // this case should never happen; if it does, there is some internal error
      throw ExceptionBlockIP(NOIJA_GENERAL, __FILE__, __LINE__);
    }
    break;
  case NETWORK:
    network_to_general();
    break;
  case IDENTITY:
    identity_to_general();
    break;
  case DIAGONAL:
    diagonal_to_general();
    break;
  case IDTY_IDTY:
    idty_idty_to_general();
    break;
  case DIAG_DIAG:
    diag_diag_to_general();
    break;
  }
  if (delete_native_format) {
    type= GENERAL;
    delete[] src;
    delete[] dst;
    delete[] d1;
    delete[] d2;
    src= NULL;
    dst= NULL;
    d1= NULL;
    d2= NULL;
  }
}

inline void MatrixBlockIP::identity_to_general()
{
  restore();
  identity_to_ija_format();
  ija_to_rowwise();
}

inline void MatrixBlockIP::diagonal_to_general()
{
  restore();
  diagonal_to_ija_format();
  ija_to_rowwise();
}

inline void MatrixBlockIP::idty_idty_to_general()
{
  restore();
  idty_idty_to_ija_format();
  ija_to_rowwise();
}

inline void MatrixBlockIP::diag_diag_to_general()
{
  restore();
  diag_diag_to_ija_format();
  ija_to_rowwise();
}

inline void MatrixBlockIP::mul_Mv_row_wise(double vout[], const double vin[])

{
  int i,j;

  for(i=0;i<m;++i){
    vout[i]=0.0;
    for(j=inirowa[i];j<=inirowa[i+1]-1;++j) vout[i]+= a[j]*vin[icola[j]];
  }
}     

inline void MatrixBlockIP::mul_Mtv_row_wise(double vout[], const double vin[])

{ 
  double r;
  int i,j;
  
  for(i=0;i<n;++i) vout[i]= 0.0;
  for(i=0;i<m;++i){
    r= vin[i];
    for(j=inirowa[i];j<=inirowa[i+1]-1;++j) vout[icola[j]]+= r*a[j];
  }
}   

inline void MatrixBlockIP::add_mul_Mv_row_wise(double vout[], const double vin[])

{
  int i,j;

  for(i=0;i<m;++i){
    for(j=inirowa[i];j<=inirowa[i+1]-1;++j) vout[i]+= a[j]*vin[icola[j]];
  }
}     

inline void MatrixBlockIP::add_mul_Mtv_row_wise(double vout[], const double vin[])

{ 
  double r;
  int i,j;
  
  for(i=0;i<m;++i){
    r= vin[i];
    for(j=inirowa[i];j<=inirowa[i+1]-1;++j) vout[icola[j]]+= r*a[j];
  }
}

inline void MatrixBlockIP::mul_Mv_column_wise(double vout[], const double vin[])

{
  double r;
  int i,j;
  
  for(i=0;i<m;++i) vout[i]= 0.0;
  for(i=0;i<n;++i){
    r= vin[i];
    for(j=inicola[i];j<inicola[i+1];++j) vout[irowa[j]]+= r*a[j];
  }
}     

inline void MatrixBlockIP::mul_Mtv_column_wise(double vout[], const double vin[])

{ 
  int i,j;

  for(i=0;i<n;++i){
    vout[i]=0.0;
    for(j=inicola[i];j<inicola[i+1];++j) vout[i]+= a[j]*vin[irowa[j]];
  }
}   

inline void MatrixBlockIP::add_mul_Mv_column_wise(double vout[], const double vin[])

{
  double r;
  int i,j;
  
  for(i=0;i<n;++i){
    r= vin[i];
    for(j=inicola[i];j<inicola[i+1];++j) vout[irowa[j]]+= r*a[j];
  }
}     

inline void MatrixBlockIP::add_mul_Mtv_column_wise(double vout[], const double vin[])

{ 
  int i,j;

  for(i=0;i<n;++i){
    for(j=inicola[i];j<inicola[i+1];++j) vout[i]+= a[j]*vin[irowa[j]];
  }
}

inline void MatrixBlockIP::mul_Mv_network(double vout[], const double vin[])

{
  double backvoutm= vout[m]; // fictitious position num_nodes backed-up

  if (type_orientation == ORIENTED) {
    // n= num_arcs, m= num_nodes-1
    int i;
    
    for(i=0;i<m+1;i++) vout[i]= 0; //this could be avoided if vout is already set to 0 by caller
    
    for(i=0;i<num_arcs;++i) {
      vout[src[i]]+= vin[i];
      vout[dst[i]]-= vin[i];
    }
  } 
  else { //NONORIENTED : matrix is [N -N]
    // n= 2*num_arcs, m= num_nodes-1
    // the following ordering of variables (arcs) is assumed:
    // first half of variables of vin[] associated to forward direction arcs;
    // second half of variables of vin[] associated to backward direction arcs
 
    int i;
    
    for(i=0;i<m+1;i++) vout[i]= 0; //this could be avoided if vout is already set to 0 by caller
    
    for(i=0;i<num_arcs;++i) {
      // forward arcs
      vout[src[i]]+= vin[i];
      vout[dst[i]]-= vin[i];
      // backward arcs
      vout[src[i]]-= vin[i+num_arcs];
      vout[dst[i]]+= vin[i+num_arcs];      
    }    
  }
  vout[m]= backvoutm; // fictitious position num_nodes restored
}

inline void MatrixBlockIP::mul_Mtv_network(double vout[], const double vin[])

{
  double * vinback;
  vinback= const_cast<double *>(vin);
  double backvinm= vin[m]; // fictitious position num_nodes backed-up

  if (type_orientation == ORIENTED) {
    // n= num_arcs, m= num_nodes-1
    int i;
    
    vinback[m]= 0.0; // fictitious position num_nodes considered
    for(i=0;i<num_arcs;++i) {
      vout[i]= 0.0;
      vout[i]+= vin[src[i]];  // fictitious position num_nodes considered
      vout[i]-= vin[dst[i]];  // fictitious position num_nodes considered
    }
  } 
  else { //NONORIENTED : matrix is [N -N]
    // n= 2*num_arcs, m= num_nodes-1
    // the following ordering of variables (arcs) is assumed:
    // first half of vout[] associated to forward direction arcs;
    // second half of vout[] associated to backward direction arcs
    int i;
    
    vinback[m]= 0.0; // fictitious position num_nodes considered
    for(i=0;i<num_arcs;++i) {
      // this is for N*vin
      vout[i]= 0.0;
      vout[i]+= vin[src[i]];  // fictitious position num_nodes considered
      vout[i]-= vin[dst[i]];  // fictitious position num_nodes considered
      // this is for (-N)*vin: just neg vout[i] previously computed
      vout[i+num_arcs]= -vout[i]; 
    }
  }
  vinback[m]= backvinm; // fictitious position num_nodes restored
}

inline void MatrixBlockIP::add_mul_Mv_network(double vout[], const double vin[])

{
  double backvoutm= vout[m]; // fictitious position num_nodes backed-up

  if (type_orientation == ORIENTED) {
    // n= num_arcs, m= num_nodes-1
    int i;
    
   // caution: fictitious position num_nodes considered in vout
    for(i=0;i<num_arcs;++i) {
      vout[src[i]]+= vin[i];
      vout[dst[i]]-= vin[i];
    }
  } 
  else { //NONORIENTED: matrix is [N -N]
    // n= 2*num_arcs, m= num_nodes-1
    // the following ordering of variables (arcs) is assumed:
    // first half of variables of vin[] associated to forward direction arcs;
    // second half of variables of vin[] associated to backward direction arcs
 
    int i;
    
   // caution: fictitious position num_nodes considered in vout
    for(i=0;i<num_arcs;++i) {
      // forward arcs
      vout[src[i]]+= vin[i];
      vout[dst[i]]-= vin[i];
      // backward arcs
      vout[src[i]]-= vin[i+num_arcs];
      vout[dst[i]]+= vin[i+num_arcs];      
    }
  }
  vout[m]= backvoutm; // fictitious position num_nodes restored
}

inline void MatrixBlockIP::add_mul_Mtv_network(double vout[], const double vin[])

{
  double * vinback;
  vinback= const_cast<double *>(vin);
  double backvinm= vin[m]; // fictitious position num_nodes backed-up

  if (type_orientation == ORIENTED) {
    // n= num_arcs, m= num_nodes-1
    int i;
    
    vinback[m]= 0.0; // fictitious position num_nodes considered
    for(i=0;i<num_arcs;++i) {
      vout[i]+= vin[src[i]];  // fictitious position num_nodes considered
      vout[i]-= vin[dst[i]];  // fictitious position num_nodes considered
    }
  } 
  else { //NONORIENTED: matrix is [N -N]
    // n= 2*num_arcs, m= num_nodes-1
    // the following ordering of variables (arcs) is assumed:
    // first half of vout[] associated to forward direction arcs;
    // second half of vout[] associated to backward direction arcs
    int i;
    
    vinback[m]= 0.0; // fictitious position num_nodes considered
    for(i=0;i<num_arcs;++i) {
      double raux= vin[src[i]]-vin[dst[i]]; // fictitious position num_nodes considered
      // this is for N*vin
      vout[i]+= raux;
      // this is for (-N)*vin
      vout[i+num_arcs]= -raux; 
    }
  }
  vinback[m]= backvinm; // fictitious position num_nodes restored
}

inline void MatrixBlockIP::mul_Mv_identity(double vout[], const double vin[])

{
  for (int i= 0; i < n; ++i) vout[i]= vin[i];
}

inline void MatrixBlockIP::add_mul_Mv_identity(double vout[], const double vin[])

{
  for (int i= 0; i < n; ++i) vout[i]+= vin[i];
}

inline void MatrixBlockIP::mul_Mv_diagonal(double vout[], const double vin[])

{
  for (int i= 0; i < n; ++i) vout[i]= d1[i]*vin[i];
}


inline void MatrixBlockIP::add_mul_Mv_diagonal(double vout[], const double vin[])

{
  for (int i= 0; i < n; ++i) vout[i]+= d1[i]*vin[i];
}


inline void MatrixBlockIP::mul_Mv_idty_idty(double vout[], const double vin[])

{
  for (int i= 0, j=m; i<m; i++, j++) vout[i]= vin[i]+vin[j];
}     

inline void MatrixBlockIP::mul_Mtv_idty_idty(double vout[], const double vin[])

{ 
  for (int i= 0, j=m; i<m; i++, j++) {
    vout[i]= vin[i];
    vout[j]= vin[i];
  }
}   

inline void MatrixBlockIP::add_mul_Mv_idty_idty(double vout[], const double vin[])

{
  for (int i= 0, j=m; i<m; i++, j++) vout[i]+= vin[i]+vin[j];
}     

inline void MatrixBlockIP::add_mul_Mtv_idty_idty(double vout[], const double vin[])

{ 
  for (int i= 0, j=m; i<m; i++, j++) {
    vout[i]+= vin[i];
    vout[j]+= vin[i];
  }
}   

inline void MatrixBlockIP::mul_Mv_diag_diag(double vout[], const double vin[])

{
  for (int i= 0, j=m; i<m; i++, j++) vout[i]= d1[i]*vin[i]+d2[i]*vin[j];
}     

inline void MatrixBlockIP::mul_Mtv_diag_diag(double vout[], const double vin[])

{ 
  for (int i= 0, j=m; i<m; i++, j++) {
    vout[i]= d1[i]*vin[i];
    vout[j]= d2[i]*vin[i];
  }
}   

inline void MatrixBlockIP::add_mul_Mv_diag_diag(double vout[], const double vin[])

{
  for (int i= 0, j=m; i<m; i++, j++) vout[i]+= d1[i]*vin[i]+d2[i]*vin[j];
}     

inline void MatrixBlockIP::add_mul_Mtv_diag_diag(double vout[], const double vin[])

{ 
  for (int i= 0, j=m; i<m; i++, j++) {
    vout[i]+= d1[i]*vin[i];
    vout[j]+= d2[i]*vin[i];
  }
}

inline void MatrixBlockIP::mul_Mv_gen_sym_uptr(double vout[], const double vin[])

{
  for(int i=0;i<m;i++) vout[i]= 0.0; 
  for(int i=0;i<m;i++){
    for(int l=inirowa[i];l<=inirowa[i+1]-1;l++) {
      int j= icola[l];
      vout[i]+= a[l]*vin[j]; // vout[i] += a[i,j]*vin[j]
      if (i != j) vout[j]+= a[l]*vin[i]; // if upper-(non-)diagonal element then vout[j] += a[j,i]*vin[i]
    }
  }
}     

inline void MatrixBlockIP::add_mul_Mv_gen_sym_uptr(double vout[], const double vin[])

{
  for(int i=0;i<m;i++){
    for(int l=inirowa[i];l<=inirowa[i+1]-1;l++) {
      int j= icola[l];
      vout[i]+= a[l]*vin[j]; // vout[i] += a[i,j]*vin[j]
      if (i != j) vout[j]+= a[l]*vin[i]; // if upper-(non-)diagonal element then vout[j] += a[j,i]*vin[i]
    }
  }
}     

inline void MatrixBlockIP::compute_D(int numBlocks, bool sameL, MatrixBlockIP L[], bool isActive[], int iniTheta[], double Theta[])

{
  if (!made_analyze_D) {
    analyze_D(numBlocks,sameL,L);
    made_analyze_D= true;
  }

  if (type==DIAGONAL) 
    compute_D_diagonal(numBlocks,sameL,L,isActive,iniTheta,Theta);
  else if (type== GEN_SYM_UPTR) 
    compute_D_gen_sym_uptr(numBlocks,sameL,L,isActive,iniTheta,Theta);
}

inline void MatrixBlockIP::analyze_D(int numBlocks, bool sameL, MatrixBlockIP L[])

{
  reset();
  m= L[0].m; // All L_i have the same number of rows
  n= m;
  
  int sizeL= sameL ? 1 : numBlocks;
  bool is_diag= true;
  for(int i=0 ; i<sizeL && is_diag ;i++) {
    if (!(L[i].type == IDENTITY || L[i].type== IDTY_IDTY || 
          L[i].type== DIAGONAL  || L[i].type== DIAG_DIAG )) is_diag= false;
  }
  if (is_diag) {
    // D is diagonal positive semidefinite matrix
    type= DIAGONAL;
    analyze_D_diagonal(numBlocks, sameL, L);
  }
  else {
    // D is symmetric positive semidefinite matrix 
    type= GEN_SYM_UPTR;
    rowwise= true;
    analyze_D_gen_sym_uptr(numBlocks, sameL, L);
  }
}


#endif //MATRIXBLOCKIP_H