BlockIPPlatform
 All Classes Files Functions Variables
MatrixBlockIP.h
Go to the documentation of this file.
1 
6 #ifndef MATRIXBLOCKIP_H
7 #define MATRIXBLOCKIP_H
8 
9 #include "SparseChol.h"
10 #include <vector>
11 #include <fstream>
12 
13 using namespace std;
14 
17  public:
18 
19  /*
20  For every type of matrix we need routines to:
21  - multiply matrix by vector, and matrix transposed by vector
22  - solve systems with matrix (either through Cholesky or specialized routines)
23  - routine to transform matrix to a GENERAL matrix (to get full matrix of constraints)
24  */
25 
27  TYPE_MATRIX type;
28  int m;
29  int n;
30  bool ija;
31  bool rowwise;
32  bool columnwise;
33 
35  int nz;
36  double *a;
37 
39  int *inirowa;
40  int *icola;
41 
43  int *inicola;
44  int *irowa;
45 
47  TYPE_ORIENTATION type_orientation;
48  int num_arcs;
49  int num_nodes;
50  int *src;
51  int *dst;
52 
54  double *d1;
55  double *d2;
56 
59 
61  int sizeL;
62  int *inirowLLt;
63  // The following vectors are allocated with ALLOC and freed with FREE
64  int *icolLLt;
65  int *blockD;
66  int *inivalD;
67  int *icolD;
68  double *valD;
69 
71  bool made_symbfct_MMt; // for symbolic factorization of M*Theta*M'
72  bool made_symbfct_M; // for symbolic factorization of M
75 
81  int num_blocks;
82 
83  // Constructor
84  MatrixBlockIP(int blocks=1);
85  // Copy constructor
87  // Destructor
88  ~MatrixBlockIP();
89 
90  // Copy
91  void copy(MatrixBlockIP *mbip);
92 
93  // Comes back to the initial state
94  void reset();
95 
96  // Comes back to the state after create the matrix
97  void restore();
98 
99  /* Native methods to create different matrix types */
100  // Native method to create a general row-wise packed matrix
101  void create_general_matrix_row_wise(int m, int n, int nz, int* &inirowa, int* &icola, double* &a);
102  // Native method to create a general column-wise packed matrix
103  void create_general_matrix_column_wise(int m, int n, int nz, int* &inicola, int* &irowa, double* &a);
104  // Native method to create a network matrix
105  void create_network_matrix(int num_arcs, int num_nodes, int* &src, int* &dst, bool oriented= true);
106  // Native method to create a identity matrix
107  void create_identity_matrix(int dim);
108  // Native method to create a identity-identity matrix
109  void create_idty_idty_matrix(int nrows);
110  // Native method to create a diagonal matrix
111  void create_diagonal_matrix(int dim, double* &d);
112  // Native method to create a diagonal-diagonal matrix
113  void create_diag_diag_matrix(int nrows, double* &d1, double* &d2);
114 /*
115  // Native method to load a general row-wise packed matrix
116  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);
117  // Native method to load a general column-wise packed matrix
118  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);
119  // Native method to load a network matrix
120  void load_network_matrix(int num_arcs, int num_nodes, int src[], int dst[], bool oriented= true, bool copy_vectors=true);
121  // Native method to load a identity matrix
122  void load_identity_matrix(int dim);
123  // Native method to load a identity-identity matrix
124  void load_idty_idty_matrix(int dim);
125  // Native method to load a diagonal matrix
126  void load_diagonal_matrix(int dim, double d[], bool copy_vectors=true);
127  // Native method to load a diagonal-diagonal matrix
128  void load_diag_diag_matrix(int dim, double d1[], double d2[], bool copy_vectors=true);
129 */
130 
131  /* Alternative methods to load/build matrices */
132  // Creates a general matrix in format ija
133  void create_general_matrix_format_ija(int m, int n, int nz, int* &row_index, int* &col_index, double* &values);
134 /*
135  // Loads a general matrix in format ija
136  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);
137 */
138  // Builds a packed rowwise format structured matrix based on diagonal blocks and linking constraints blocks
139  void compute_full_matrix(int numBlocks, bool sameN, MatrixBlockIP N[], bool sameL, MatrixBlockIP L[], int numActiveLnk=0, int *listActiveLnk=NULL);
140  // Create internal structure arrays to compute D = Theta_(numBlocks+1) + sum{i..numBlocks} L_i*Theta_i*L_i'
141  void analyze_D(int numBlocks, bool sameL, MatrixBlockIP L[]);
142  void analyze_D_diagonal(int numBlocks, bool sameL, MatrixBlockIP L[]);
143  void analyze_D_gen_sym_uptr(int numBlocks, bool sameL, MatrixBlockIP L[]);
144  // Compute D = Theta_(numBlocks+1) + sum{i..numBlocks} L_i*Theta_i*L_i'
145  void compute_D(int numBlocks, bool sameL, MatrixBlockIP L[], bool isActive[], int iniTheta[], double Theta[]);
146  void compute_D_diagonal(int numBlocks, bool sameL, MatrixBlockIP L[], bool isActive[], int iniTheta[], double Theta[]);
147  void compute_D_gen_sym_uptr(int numBlocks, bool sameL, MatrixBlockIP L[], bool isActive[], int iniTheta[], double Theta[]);
148 
149 
150  /* Matrix transformations */
151  // Calculates and stores the matrix in packed rowwise format
152  void native_to_general(bool delete_native_format=false);
153  // Calculates and stores the matrix in packed rowwise format
154  void ija_to_rowwise();
155  // Calculates and stores the network matrix (either oriented or nonoriented) in packed rowwise format
156  void network_to_general();
157  // Calculates and stores the I matrix in packed rowwise format
158  void identity_to_general();
159  // Calculates and stores the diagonal D matrix in packed rowwise format
160  void diagonal_to_general();
161  // Calculates and stores the matrix [I I] in packed rowwise format
162  void idty_idty_to_general();
163  // Calculates and stores the matrix [D1 D2] in packed rowwise format
164  void diag_diag_to_general();
165 
166  // Calculates and stores the network matrix in ija format.
167  void network_to_ija_format();
168  // Calculates and stores the identity I matrix in ija format.
169  void identity_to_ija_format();
170  // Calculates and stores the diagonal D1 matrix in ija format.
171  void diagonal_to_ija_format();
172  // Calculates and stores the matrix [I I] in ija format
173  void idty_idty_to_ija_format();
174  // Calculates and stores the matrix [D1 D2] in ija format
175  void diag_diag_to_ija_format();
176  // Order the matrix when has been created in row-wise or column-wise format
177  void order_matrix();
178 
179  /* Matrix operations */
180  // Matrix-vector product vout=M*vin (vout(m),vin(n))
181  void mul_Mv(double vout[], const double vin[]);
182  // MatrixTranspose-vector product vout=M(t)*vin (vout(n),vin(m)
183  void mul_Mtv(double vout[], const double vin[]);
184  // Add matrix-vector product vout += M*vin (vout(m),vin(n))
185  void add_mul_Mv(double vout[], const double vin[]);
186  // Add matrixTranspose-vector product vout += M(t)*vin (vout(n),vin(m)
187  void add_mul_Mtv(double vout[], const double vin[]);
188  // Given a column determines if exists an element for each row of the matrix, in case of exist returns the value.
189  void exist_var_in_row(int column, bool appear[], double coefs[]);
190  // Adds a new column into the matrix
191  void add_new_column(int size, int irows[], double values[]);
192  // Adds new columns into the matrix
193  void add_new_columns(int num_columns, int size[], int *irows[], double *values[]);
194  // Changes the sign of some columns
195  void change_columns_sign(int size, int columns[]);
196  // Changes the sign of some rows
197  void change_rows_sign(int size, int rows[]);
198  // Deletes some rows
199  void delete_rows(int size, int rows[]);
200  // Gets a column of the matrix
201  int get_column(int column, int* &irows, double* &values, bool invert_sign= false);
202 
203  // SparseChol interface functions
204  // next three for systems with M*Theta*M'
205  void symbolic_fact_MMt(CHOL_SOLVER chslv=SPRSBLKLLT, int *prov_pfa= NULL, int *prov_ipfa= NULL, NUMBERING prov_pfa_numbering= NOT_COMPUTED);
206  void numeric_fact_MMt(double *Theta, int i_k=0);
207  void numeric_solve_MMt(double *rhs, int i_k= 0, WHO_PERMUTES whoperm= CHOLESKY);
208  // next three for systems with M (M positive semidefinite matrix)
209  void symbolic_fact_M(CHOL_SOLVER chslv=SPRSBLKLLT, int *prov_pfa= NULL, int *prov_ipfa= NULL, NUMBERING prov_pfa_numbering= NOT_COMPUTED);
210  void numeric_fact_M();
211  void numeric_solve_M(double *rhs, WHO_PERMUTES whoperm= CHOLESKY);
212 
213  // access to data in class
214  int get_pfa(int i); // returns SparseChol pfa[i] with no check
215  int get_ipfa(int i); // returns SparseChol ipfa[i] with no check
216  int get_maxlnz(); // returns SparseChol maxlnz with no check
217  int get_maxfillin(); // returns SparseChol maxfillin with no check
218  int get_njka(); // returns SparseChol njka with no check
219  int get_num_zero_pivots(); // returns SparseChol num_zero_pivots with no check
220  int get_num_semidef_matrix(); // returns SparseChol num_semidef_matrix with no check
221 
222  // only to test
223  void print_matrix(bool ija= true, bool start_one= true);
224  void print_matrix(ofstream &outfile, bool print_ija= true, bool start_one= true);
225  // Print some positions and name of a vector
226  void print_vector(double *v, int size, string name);
227 
228  // Pre: icola is ordered
229  // Convert a matrix in column-wise format to row-wise format.
230  void column_wise_to_row_wise_format();
231 
232  // Pre: irowa is ordered
233  // Convert a matrix in row-wise format to column-wise format.
234  void row_wise_to_column_wise_format();
235 
236  private:
237 
238  // free all the memory allocated by the application - not the user
239  void free_memory();
240 
241  void mul_Mv_row_wise(double vout[], const double vin[]);
242  void mul_Mtv_row_wise(double vout[], const double vin[]);
243  void add_mul_Mv_row_wise(double vout[], const double vin[]);
244  void add_mul_Mtv_row_wise(double vout[], const double vin[]);
245 
246  void mul_Mv_column_wise(double vout[], const double vin[]);
247  void mul_Mtv_column_wise(double vout[], const double vin[]);
248  void add_mul_Mv_column_wise(double vout[], const double vin[]);
249  void add_mul_Mtv_column_wise(double vout[], const double vin[]);
250 
251  void mul_Mv_network(double vout[], const double vin[]);
252  void mul_Mtv_network(double vout[], const double vin[]);
253  void add_mul_Mv_network(double vout[], const double vin[]);
254  void add_mul_Mtv_network(double vout[], const double vin[]);
255 
256  void mul_Mv_identity(double vout[], const double vin[]);
257  void add_mul_Mv_identity(double vout[], const double vin[]);
258 
259  void mul_Mv_diagonal(double vout[], const double vin[]);
260  void add_mul_Mv_diagonal(double vout[], const double vin[]);
261 
262  void mul_Mv_idty_idty(double vout[], const double vin[]);
263  void mul_Mtv_idty_idty(double vout[], const double vin[]);
264  void add_mul_Mv_idty_idty(double vout[], const double vin[]);
265  void add_mul_Mtv_idty_idty(double vout[], const double vin[]);
266 
267  void mul_Mv_diag_diag(double vout[], const double vin[]);
268  void mul_Mtv_diag_diag(double vout[], const double vin[]);
269  void add_mul_Mv_diag_diag(double vout[], const double vin[]);
270  void add_mul_Mtv_diag_diag(double vout[], const double vin[]);
271 
272  void mul_Mv_gen_sym_uptr(double vout[], const double vin[]); // M=M', then mul_Mtv call this one
273  void add_mul_Mv_gen_sym_uptr(double vout[], const double vin[]); // M=M', then mul_Mtv call this one
274 
275  // Order a matrix packed format (both column-wise and row-wise)
276  void order_packed_format(int begsize, int nz, int* &beg, int* &ind, double* &val);
277 
279  struct Order_ija {
280  int *row;
281  int *col;
282  bool operator() (int i,int j)
283  {
284  return ((row[i] < row[j]) || ((row[i] == row[j]) && (col[i] < col[j])));
285  }
286  };
287 
289  struct Order_vector {
290  int *v;
291  bool operator() (int i,int j)
292  {
293  return (v[i] < v[j]);
294  }
295  };
296 };
297 
298 
299 inline void MatrixBlockIP::symbolic_fact_MMt(CHOL_SOLVER chslv, int *prov_pfa, int *prov_ipfa, NUMBERING prov_pfa_numbering)
301 {
302  chol.initialize(chslv,type,type_orientation);
303 
304  if (m <= 0) return; // for 0 rows matrices do nothing
305 
306  switch (type) {
307  case GENERAL:
308  chol.symbolic_fact_MMt(m, n, nz, icola, inirowa, a, num_blocks, prov_pfa, prov_ipfa, prov_pfa_numbering);
309  break;
310  case NETWORK:
311  chol.symbolic_fact_MMt(num_nodes, num_arcs, src, dst, num_blocks, prov_pfa, prov_ipfa, prov_pfa_numbering);
312  break;
313  case IDENTITY:
314  case IDTY_IDTY:
315  chol.symbolic_fact_MMt(m, num_blocks);
316  break;
317  case DIAGONAL:
318  chol.symbolic_fact_MMt(m,d1,num_blocks);
319  break;
320  case DIAG_DIAG:
321  chol.symbolic_fact_MMt(m,d1,d2,num_blocks);
322  break;
323  default:
324  throw ExceptionBlockIP(NOT_FOR_THIS_MATRIX_TYPE, __FILE__, __LINE__);
325  }
326 }
327 
328 inline void MatrixBlockIP::numeric_fact_MMt(double *Theta, int i_k)
330 {
331  if (m <= 0) return; // for 0 rows matrices do nothing
332 
333  // if first time then creation and symbolic factorization of M*Theta*M'
334  if (!made_symbfct_MMt) {
335  symbolic_fact_MMt();
336  made_symbfct_MMt= true;
337  }
338 
339  chol.numeric_fact_MMt(Theta, i_k);
340 }
341 
342 inline void MatrixBlockIP::numeric_solve_MMt(double *rhs, int i_k, WHO_PERMUTES whoperm)
344 {
345  if (m <= 0) return; // for 0 rows matrices do nothing
346 
347  chol.numeric_solve_MMt(rhs, i_k, whoperm);
348 }
349 
350 inline void MatrixBlockIP::symbolic_fact_M(CHOL_SOLVER chslv, int *prov_pfa, int *prov_ipfa, NUMBERING prov_pfa_numbering)
352 {
353  chol.initialize(chslv,type,type_orientation);
354 
355  if (m <= 0) return; // for 0 rows matrices do nothing
356 
357  // only implemented (and needed by now) for GEN_SYM_UPTR and DIAGONAL matrices
358  switch (type) {
359  case DIAGONAL:
360  chol.symbolic_fact_M(m);
361  break;
362  case GEN_SYM_UPTR:
363  chol.symbolic_fact_M(m, nz, icola, inirowa, prov_pfa, prov_ipfa, prov_pfa_numbering);
364  break;
365  default:
366  throw ExceptionBlockIP(NOT_FOR_THIS_MATRIX_TYPE, __FILE__, __LINE__);
367  }
368 }
369 
372 {
373  if (m <= 0) return; // for 0 rows matrices do nothing
374 
375  // if first time then symbolic factorization of M
376  if (!made_symbfct_M) {
377  symbolic_fact_M();
378  made_symbfct_M= true;
379  }
380 
381  switch (type) {
382  case DIAGONAL:
383  chol.numeric_fact_M(d1);
384  break;
385  case GEN_SYM_UPTR:
386  chol.numeric_fact_M(a);
387  break;
388  default:
389  throw ExceptionBlockIP(NOT_FOR_THIS_MATRIX_TYPE, __FILE__, __LINE__);
390  }
391 }
392 
393 inline void MatrixBlockIP::numeric_solve_M(double *rhs, WHO_PERMUTES whoperm)
395 {
396  if (m <= 0) return; // for 0 rows matrices do nothing
397 
398  chol.numeric_solve_M(rhs, whoperm);
399 }
400 
401 inline int MatrixBlockIP::get_pfa(int i)
402 
403 
407 {
408  return(chol.get_pfa(i));
409 }
410 
411 inline int MatrixBlockIP::get_ipfa(int i)
413 
417 {
418  return(chol.get_ipfa(i));
419 }
420 
423 
427 {
428  return (chol.get_maxlnz());
429 }
430 
433 
437 {
438  return (chol.get_maxfillin());
439 }
440 
443 
447 {
448  return (chol.get_njka());
449 }
450 
453 
457 {
458  return (chol.get_num_zero_pivots());
459 }
460 
463 
467 {
468  return (chol.get_num_semidef_matrix());
469 }
470 
471 inline void MatrixBlockIP::native_to_general(bool delete_native_format)
473 
476 {
477  if (rowwise) return; // already stored
478 
479  switch (type) {
480  case GENERAL:
481  if (columnwise)
482  column_wise_to_row_wise_format();
483  else if (ija)
484  ija_to_rowwise();
485  else {
486  // this case should never happen; if it does, there is some internal error
487  throw ExceptionBlockIP(NOIJA_GENERAL, __FILE__, __LINE__);
488  }
489  break;
490  case NETWORK:
491  network_to_general();
492  break;
493  case IDENTITY:
494  identity_to_general();
495  break;
496  case DIAGONAL:
497  diagonal_to_general();
498  break;
499  case IDTY_IDTY:
500  idty_idty_to_general();
501  break;
502  case DIAG_DIAG:
503  diag_diag_to_general();
504  break;
505  }
506  if (delete_native_format) {
507  type= GENERAL;
508  delete[] src;
509  delete[] dst;
510  delete[] d1;
511  delete[] d2;
512  src= NULL;
513  dst= NULL;
514  d1= NULL;
515  d2= NULL;
516  }
517 }
518 
521 {
522  restore();
523  identity_to_ija_format();
524  ija_to_rowwise();
525 }
526 
529 {
530  restore();
531  diagonal_to_ija_format();
532  ija_to_rowwise();
533 }
534 
537 {
538  restore();
539  idty_idty_to_ija_format();
540  ija_to_rowwise();
541 }
542 
545 {
546  restore();
547  diag_diag_to_ija_format();
548  ija_to_rowwise();
549 }
550 
551 inline void MatrixBlockIP::mul_Mv_row_wise(double vout[], const double vin[])
553 
557 {
558  int i,j;
559 
560  for(i=0;i<m;++i){
561  vout[i]=0.0;
562  for(j=inirowa[i];j<=inirowa[i+1]-1;++j) vout[i]+= a[j]*vin[icola[j]];
563  }
564 }
565 
566 inline void MatrixBlockIP::mul_Mtv_row_wise(double vout[], const double vin[])
568 
572 {
573  double r;
574  int i,j;
575 
576  for(i=0;i<n;++i) vout[i]= 0.0;
577  for(i=0;i<m;++i){
578  r= vin[i];
579  for(j=inirowa[i];j<=inirowa[i+1]-1;++j) vout[icola[j]]+= r*a[j];
580  }
581 }
582 
583 inline void MatrixBlockIP::add_mul_Mv_row_wise(double vout[], const double vin[])
585 
589 {
590  int i,j;
591 
592  for(i=0;i<m;++i){
593  for(j=inirowa[i];j<=inirowa[i+1]-1;++j) vout[i]+= a[j]*vin[icola[j]];
594  }
595 }
596 
597 inline void MatrixBlockIP::add_mul_Mtv_row_wise(double vout[], const double vin[])
599 
603 {
604  double r;
605  int i,j;
606 
607  for(i=0;i<m;++i){
608  r= vin[i];
609  for(j=inirowa[i];j<=inirowa[i+1]-1;++j) vout[icola[j]]+= r*a[j];
610  }
611 }
612 
613 inline void MatrixBlockIP::mul_Mv_column_wise(double vout[], const double vin[])
615 
619 {
620  double r;
621  int i,j;
622 
623  for(i=0;i<m;++i) vout[i]= 0.0;
624  for(i=0;i<n;++i){
625  r= vin[i];
626  for(j=inicola[i];j<inicola[i+1];++j) vout[irowa[j]]+= r*a[j];
627  }
628 }
629 
630 inline void MatrixBlockIP::mul_Mtv_column_wise(double vout[], const double vin[])
632 
636 {
637  int i,j;
638 
639  for(i=0;i<n;++i){
640  vout[i]=0.0;
641  for(j=inicola[i];j<inicola[i+1];++j) vout[i]+= a[j]*vin[irowa[j]];
642  }
643 }
644 
645 inline void MatrixBlockIP::add_mul_Mv_column_wise(double vout[], const double vin[])
647 
651 {
652  double r;
653  int i,j;
654 
655  for(i=0;i<n;++i){
656  r= vin[i];
657  for(j=inicola[i];j<inicola[i+1];++j) vout[irowa[j]]+= r*a[j];
658  }
659 }
660 
661 inline void MatrixBlockIP::add_mul_Mtv_column_wise(double vout[], const double vin[])
663 
667 {
668  int i,j;
669 
670  for(i=0;i<n;++i){
671  for(j=inicola[i];j<inicola[i+1];++j) vout[i]+= a[j]*vin[irowa[j]];
672  }
673 }
674 
675 inline void MatrixBlockIP::mul_Mv_network(double vout[], const double vin[])
677 
684 {
685  double backvoutm= vout[m]; // fictitious position num_nodes backed-up
686 
687  if (type_orientation == ORIENTED) {
688  // n= num_arcs, m= num_nodes-1
689  int i;
690 
691  for(i=0;i<m+1;i++) vout[i]= 0; //this could be avoided if vout is already set to 0 by caller
692 
693  for(i=0;i<num_arcs;++i) {
694  vout[src[i]]+= vin[i];
695  vout[dst[i]]-= vin[i];
696  }
697  }
698  else { //NONORIENTED : matrix is [N -N]
699  // n= 2*num_arcs, m= num_nodes-1
700  // the following ordering of variables (arcs) is assumed:
701  // first half of variables of vin[] associated to forward direction arcs;
702  // second half of variables of vin[] associated to backward direction arcs
703 
704  int i;
705 
706  for(i=0;i<m+1;i++) vout[i]= 0; //this could be avoided if vout is already set to 0 by caller
707 
708  for(i=0;i<num_arcs;++i) {
709  // forward arcs
710  vout[src[i]]+= vin[i];
711  vout[dst[i]]-= vin[i];
712  // backward arcs
713  vout[src[i]]-= vin[i+num_arcs];
714  vout[dst[i]]+= vin[i+num_arcs];
715  }
716  }
717  vout[m]= backvoutm; // fictitious position num_nodes restored
718 }
719 
720 inline void MatrixBlockIP::mul_Mtv_network(double vout[], const double vin[])
722 
729 {
730  double * vinback;
731  vinback= const_cast<double *>(vin);
732  double backvinm= vin[m]; // fictitious position num_nodes backed-up
733 
734  if (type_orientation == ORIENTED) {
735  // n= num_arcs, m= num_nodes-1
736  int i;
737 
738  vinback[m]= 0.0; // fictitious position num_nodes considered
739  for(i=0;i<num_arcs;++i) {
740  vout[i]= 0.0;
741  vout[i]+= vin[src[i]]; // fictitious position num_nodes considered
742  vout[i]-= vin[dst[i]]; // fictitious position num_nodes considered
743  }
744  }
745  else { //NONORIENTED : matrix is [N -N]
746  // n= 2*num_arcs, m= num_nodes-1
747  // the following ordering of variables (arcs) is assumed:
748  // first half of vout[] associated to forward direction arcs;
749  // second half of vout[] associated to backward direction arcs
750  int i;
751 
752  vinback[m]= 0.0; // fictitious position num_nodes considered
753  for(i=0;i<num_arcs;++i) {
754  // this is for N*vin
755  vout[i]= 0.0;
756  vout[i]+= vin[src[i]]; // fictitious position num_nodes considered
757  vout[i]-= vin[dst[i]]; // fictitious position num_nodes considered
758  // this is for (-N)*vin: just neg vout[i] previously computed
759  vout[i+num_arcs]= -vout[i];
760  }
761  }
762  vinback[m]= backvinm; // fictitious position num_nodes restored
763 }
764 
765 inline void MatrixBlockIP::add_mul_Mv_network(double vout[], const double vin[])
767 
774 {
775  double backvoutm= vout[m]; // fictitious position num_nodes backed-up
776 
777  if (type_orientation == ORIENTED) {
778  // n= num_arcs, m= num_nodes-1
779  int i;
780 
781  // caution: fictitious position num_nodes considered in vout
782  for(i=0;i<num_arcs;++i) {
783  vout[src[i]]+= vin[i];
784  vout[dst[i]]-= vin[i];
785  }
786  }
787  else { //NONORIENTED: matrix is [N -N]
788  // n= 2*num_arcs, m= num_nodes-1
789  // the following ordering of variables (arcs) is assumed:
790  // first half of variables of vin[] associated to forward direction arcs;
791  // second half of variables of vin[] associated to backward direction arcs
792 
793  int i;
794 
795  // caution: fictitious position num_nodes considered in vout
796  for(i=0;i<num_arcs;++i) {
797  // forward arcs
798  vout[src[i]]+= vin[i];
799  vout[dst[i]]-= vin[i];
800  // backward arcs
801  vout[src[i]]-= vin[i+num_arcs];
802  vout[dst[i]]+= vin[i+num_arcs];
803  }
804  }
805  vout[m]= backvoutm; // fictitious position num_nodes restored
806 }
807 
808 inline void MatrixBlockIP::add_mul_Mtv_network(double vout[], const double vin[])
810 
817 {
818  double * vinback;
819  vinback= const_cast<double *>(vin);
820  double backvinm= vin[m]; // fictitious position num_nodes backed-up
821 
822  if (type_orientation == ORIENTED) {
823  // n= num_arcs, m= num_nodes-1
824  int i;
825 
826  vinback[m]= 0.0; // fictitious position num_nodes considered
827  for(i=0;i<num_arcs;++i) {
828  vout[i]+= vin[src[i]]; // fictitious position num_nodes considered
829  vout[i]-= vin[dst[i]]; // fictitious position num_nodes considered
830  }
831  }
832  else { //NONORIENTED: matrix is [N -N]
833  // n= 2*num_arcs, m= num_nodes-1
834  // the following ordering of variables (arcs) is assumed:
835  // first half of vout[] associated to forward direction arcs;
836  // second half of vout[] associated to backward direction arcs
837  int i;
838 
839  vinback[m]= 0.0; // fictitious position num_nodes considered
840  for(i=0;i<num_arcs;++i) {
841  double raux= vin[src[i]]-vin[dst[i]]; // fictitious position num_nodes considered
842  // this is for N*vin
843  vout[i]+= raux;
844  // this is for (-N)*vin
845  vout[i+num_arcs]= -raux;
846  }
847  }
848  vinback[m]= backvinm; // fictitious position num_nodes restored
849 }
850 
851 inline void MatrixBlockIP::mul_Mv_identity(double vout[], const double vin[])
853 
858 {
859  for (int i= 0; i < n; ++i) vout[i]= vin[i];
860 }
861 
862 inline void MatrixBlockIP::add_mul_Mv_identity(double vout[], const double vin[])
864 
869 {
870  for (int i= 0; i < n; ++i) vout[i]+= vin[i];
871 }
872 
873 inline void MatrixBlockIP::mul_Mv_diagonal(double vout[], const double vin[])
875 
879 {
880  for (int i= 0; i < n; ++i) vout[i]= d1[i]*vin[i];
881 }
882 
883 
884 inline void MatrixBlockIP::add_mul_Mv_diagonal(double vout[], const double vin[])
886 
890 {
891  for (int i= 0; i < n; ++i) vout[i]+= d1[i]*vin[i];
892 }
893 
894 
895 inline void MatrixBlockIP::mul_Mv_idty_idty(double vout[], const double vin[])
897 
901 {
902  for (int i= 0, j=m; i<m; i++, j++) vout[i]= vin[i]+vin[j];
903 }
904 
905 inline void MatrixBlockIP::mul_Mtv_idty_idty(double vout[], const double vin[])
907 
911 {
912  for (int i= 0, j=m; i<m; i++, j++) {
913  vout[i]= vin[i];
914  vout[j]= vin[i];
915  }
916 }
917 
918 inline void MatrixBlockIP::add_mul_Mv_idty_idty(double vout[], const double vin[])
920 
924 {
925  for (int i= 0, j=m; i<m; i++, j++) vout[i]+= vin[i]+vin[j];
926 }
927 
928 inline void MatrixBlockIP::add_mul_Mtv_idty_idty(double vout[], const double vin[])
930 
934 {
935  for (int i= 0, j=m; i<m; i++, j++) {
936  vout[i]+= vin[i];
937  vout[j]+= vin[i];
938  }
939 }
940 
941 inline void MatrixBlockIP::mul_Mv_diag_diag(double vout[], const double vin[])
943 
947 {
948  for (int i= 0, j=m; i<m; i++, j++) vout[i]= d1[i]*vin[i]+d2[i]*vin[j];
949 }
950 
951 inline void MatrixBlockIP::mul_Mtv_diag_diag(double vout[], const double vin[])
953 
957 {
958  for (int i= 0, j=m; i<m; i++, j++) {
959  vout[i]= d1[i]*vin[i];
960  vout[j]= d2[i]*vin[i];
961  }
962 }
963 
964 inline void MatrixBlockIP::add_mul_Mv_diag_diag(double vout[], const double vin[])
966 
970 {
971  for (int i= 0, j=m; i<m; i++, j++) vout[i]+= d1[i]*vin[i]+d2[i]*vin[j];
972 }
973 
974 inline void MatrixBlockIP::add_mul_Mtv_diag_diag(double vout[], const double vin[])
976 
980 {
981  for (int i= 0, j=m; i<m; i++, j++) {
982  vout[i]+= d1[i]*vin[i];
983  vout[j]+= d2[i]*vin[i];
984  }
985 }
986 
987 inline void MatrixBlockIP::mul_Mv_gen_sym_uptr(double vout[], const double vin[])
989 
993 {
994  for(int i=0;i<m;i++) vout[i]= 0.0;
995  for(int i=0;i<m;i++){
996  for(int l=inirowa[i];l<=inirowa[i+1]-1;l++) {
997  int j= icola[l];
998  vout[i]+= a[l]*vin[j]; // vout[i] += a[i,j]*vin[j]
999  if (i != j) vout[j]+= a[l]*vin[i]; // if upper-(non-)diagonal element then vout[j] += a[j,i]*vin[i]
1000  }
1001  }
1002 }
1003 
1004 inline void MatrixBlockIP::add_mul_Mv_gen_sym_uptr(double vout[], const double vin[])
1006 
1010 {
1011  for(int i=0;i<m;i++){
1012  for(int l=inirowa[i];l<=inirowa[i+1]-1;l++) {
1013  int j= icola[l];
1014  vout[i]+= a[l]*vin[j]; // vout[i] += a[i,j]*vin[j]
1015  if (i != j) vout[j]+= a[l]*vin[i]; // if upper-(non-)diagonal element then vout[j] += a[j,i]*vin[i]
1016  }
1017  }
1018 }
1019 
1020 inline void MatrixBlockIP::compute_D(int numBlocks, bool sameL, MatrixBlockIP L[], bool isActive[], int iniTheta[], double Theta[])
1022 
1029 {
1030  if (!made_analyze_D) {
1031  analyze_D(numBlocks,sameL,L);
1032  made_analyze_D= true;
1033  }
1034 
1035  if (type==DIAGONAL)
1036  compute_D_diagonal(numBlocks,sameL,L,isActive,iniTheta,Theta);
1037  else if (type== GEN_SYM_UPTR)
1038  compute_D_gen_sym_uptr(numBlocks,sameL,L,isActive,iniTheta,Theta);
1039 }
1040 
1041 inline void MatrixBlockIP::analyze_D(int numBlocks, bool sameL, MatrixBlockIP L[])
1044 
1049 {
1050  reset();
1051  m= L[0].m; // All L_i have the same number of rows
1052  n= m;
1053 
1054  int sizeL= sameL ? 1 : numBlocks;
1055  bool is_diag= true;
1056  for(int i=0 ; i<sizeL && is_diag ;i++) {
1057  if (!(L[i].type == IDENTITY || L[i].type== IDTY_IDTY ||
1058  L[i].type== DIAGONAL || L[i].type== DIAG_DIAG )) is_diag= false;
1059  }
1060  if (is_diag) {
1061  // D is diagonal positive semidefinite matrix
1062  type= DIAGONAL;
1063  analyze_D_diagonal(numBlocks, sameL, L);
1064  }
1065  else {
1066  // D is symmetric positive semidefinite matrix
1067  type= GEN_SYM_UPTR;
1068  rowwise= true;
1069  analyze_D_gen_sym_uptr(numBlocks, sameL, L);
1070  }
1071 }
1072 
1073 
1074 #endif //MATRIXBLOCKIP_H