Main Page   Namespace List   Class Hierarchy   Alphabetical List   Compound List   File List   Namespace Members   Compound Members   File Members  

matvec::SparseMatrix Class Reference

#include <sparsematrix.h>

List of all members.

Detailed Description

a sparse matrix is a matrix with very sparse elements

See also:
Matrix

Definition at line 54 of file sparsematrix.h.

Public Methods

 SparseMatrix (void)
 SparseMatrix (int n, int max_nz)
 SparseMatrix (const SparseMatrix &A)
 ~SparseMatrix (void)
const SparseMatrix & operator= (const SparseMatrix &A)
double operator() (const int i, const int j)
SparseMatrix & resize (const int n, const int maxnz)
double q (const Vector< double > &v1, const Vector< double > &v2)
double q (const double *x, const double *y)
double q (const double *x, const double *y, const int i1, const int i2)
double qTLW (const Vector< double > &v1, const Vector< double > &v2)
double qTLW (const double *x, const double *y)
double qTLW (const double *x, const double *y, const int i1, const int i2)
unsigned num_rows (void) const
unsigned nz (void) const
unsigned close (void)
unsigned factorization (const int path=5, const double stopval=1.0e-10)
double * a (void)
doubleMatrix dense (unsigned r1, unsigned r2, unsigned c1, unsigned c2) const
doubleMatrix dense (void)
Vector< double > row (const int ithrow)
double logdet (unsigned &irank)
double logdet (unsigned *irank)
double logdet (void)
unsigned irank (void)
int solve (double *x, const double *b, const std::string &method="ysmp", double relax=1.0, double stopval=0.001, int mxiter=1000)
int solve (Vector< double > &x, const Vector< double > &b, const std::string &method="ysmp", double relax=1.0, double stopval=0.001, int mxiter=1000)
int * ia (void)
int * ja (void)
int reset (const int i, const int j, const double value)
int insert (const int i, const int j, const double value)
void reorder (const int path=2)
void add_diag (const int ibeg, const int iend, const double r)
void gibbs_iterate (double *x, const double *r, double *w, const double se=1.0)
void display (unsigned r1, unsigned r2, unsigned c1, unsigned c2, const std::string &msg="")
void display (const std::string &msg="")
void release_nodestuff (void)
void save (const std::string &fname, const int io_mode=std::ios::out) const
void release (void)
void mv (const Vector< double > &v, Vector< double > &result)
void mv (const double *v, const unsigned n, double *result)
double getaij (const int i, const int j)
int initialize_node_list (void)
int minfilsymelm_node_list (void)
int elim_node (unsigned i, Vector< unsigned > Aorder)
void display_node_list (void)
void display_order (void)
void display_inv_order (void)
int initial_lij (void)
void display_lij (void)
int factor (void)
int solv (double *x, const double *rhs)
int solvrow (double *x, const double *rhs, unsigned row, Vector< double > &sol, Vector< double > &temp_rhs, unsigned *row_pt, unsigned *col_pt)
Vector< double > solv (Vector< double > rhs)

Public Attributes

double logdeterm
Vector< double > info_vec

Protected Methods

void copyfrom (const SparseMatrix &A)

Protected Attributes

Vector< NodeStruct > node_list
Vector< unsigned > mask
Vector< unsigned > masklist
Vector< unsigned > order
Vector< unsigned > inv_order
Vector< double > diag
unsigned dim
unsigned nonzero
unsigned max_nz
unsigned hsize
unsigned insertend
unsigned nsp
int * iiv
int * jjv
int * perm
int * iperm
idxtype * xadj
idxtype * adjncy
double * aav
double * rsp
int * srt_hash
int * hash_srt
char solver_name
unsigned factor_done
unsigned elim_done
unsigned initial_lij_done
unsigned initialize_node_list_done
unsigned logdet_done
unsigned rank

Friends

class GLMM
int operator== (const SparseMatrix &A, const SparseMatrix &B)
int operator!= (const SparseMatrix &A, const SparseMatrix &B)
SparseMatrix operator * (const SparseMatrix &A, const double s)
SparseMatrix operator * (const double s, const SparseMatrix &A)
Vector< double > operator * (const SparseMatrix &A, const Vector< double > &V)
std::ostream & operator<< (std::ostream &stream, const SparseMatrix &A)

Constructor & Destructor Documentation

matvec::SparseMatrix::SparseMatrix (  void    ) 
 

A constructor to create an SparseMatrix object.

Definition at line 140 of file sparsematrix.cpp.

References aav, adjncy, dim, elim_done, factor_done, hash_srt, hsize, iiv, info_vec, initial_lij_done, initialize_node_list_done, insertend, jjv, logdet_done, max_nz, nonzero, matvec::Vector< double >::reserve(), solver_name, srt_hash, and xadj. 

00141 {
00142    dim         = 0;
00143    nonzero     = 0;
00144    max_nz      = 0;
00145    hsize       = 0;
00146    insertend   = 0;
00147    factor_done = 0;
00148    elim_done   = 0;
00149    logdet_done = 0;
00150    initial_lij_done = 0;
00151    initialize_node_list_done = 0;
00152    iiv         = 0;
00153    jjv         = 0;
00154    xadj        = 0;
00155    adjncy      = 0;
00156    aav         = 0;
00157    srt_hash    = 0;
00158    hash_srt    = 0;
00159 
00160    solver_name = 'N';
00161    info_vec.reserve(2);
00162 }

matvec::SparseMatrix::SparseMatrix (  int    n,
int    maxnz
) 
 

A constructor to create an SparseMatrix object of size n by n and maximum non-zero elemenat maxnz.

Definition at line 195 of file sparsematrix.cpp.

References aav, dim, elim_done, factor_done, hash_srt, hsize, iiv, info_vec, initial_lij_done, initialize_node_list_done, insert(), insertend, jjv, logdet_done, max_nz, matvec::next_prime(), nonzero, matvec::Vector< double >::reserve(), solver_name, and srt_hash. 

00196 {
00197    if (n<0 || maxnz<0) throw exception("SparseMatrix(m,n): requring nonnegative m,n");
00198    dim = n;
00199    max_nz = maxnz;
00200    nonzero = 0;
00201    factor_done = 0;
00202    elim_done = 0;
00203    logdet_done = 0;
00204    initial_lij_done = 0;
00205    initialize_node_list_done = 0;
00206    if (maxnz > 0) {
00207       hsize   = static_cast<unsigned>(next_prime(static_cast<long>(1.25*maxnz + 50)));
00208    }
00209    iiv     = new int [hsize];
00210    jjv     = new int [hsize];
00211    aav     = new double [hsize];
00212    hash_srt = new int [hsize+1];
00213    srt_hash = new int [hsize+1];
00214    nonzero=0;
00215    solver_name  = 'N';
00216    insertend = 0;
00217    info_vec.reserve(2);
00218 
00219    memset(iiv,'\0',sizeof(int)*hsize);
00220    memset(jjv,'\0',sizeof(int)*hsize);
00221    memset(aav,'\0',sizeof(double)*hsize);
00222    if(hsize){
00223      memset(hash_srt,'\0',sizeof(int)*(hsize+1));
00224      memset(srt_hash,'\0',sizeof(int)*(hsize+1));
00225    }
00226    nonzero=0;
00227    for (unsigned i=1; i<=dim; i++) insert(i,i,0.0);
00228 }

matvec::SparseMatrix::SparseMatrix (  const SparseMatrix &    A ) 
 

A constructor to create an identical SparseMatrix object to A.

Definition at line 167 of file sparsematrix.cpp.

References aav, adjncy, copyfrom(), dim, elim_done, factor_done, hash_srt, hsize, iiv, initial_lij_done, initialize_node_list_done, insertend, jjv, logdet_done, max_nz, nonzero, solver_name, srt_hash, and xadj. 

00168 {
00169    dim         = 0;
00170    nonzero     = 0;
00171    max_nz      = 0;
00172    hsize       = 0;
00173    insertend   = 0;
00174    factor_done = 0;
00175    elim_done   = 0;
00176    logdet_done = 0;
00177    initial_lij_done = 0;  
00178    initialize_node_list_done = 0;
00179    iiv         = 0;
00180    jjv         = 0;
00181    xadj        = 0;
00182    adjncy      = 0;
00183    aav         = 0;
00184    srt_hash    = 0;
00185    hash_srt    = 0;
00186 
00187    solver_name  = 'N';
00188    copyfrom(A);
00189 }

matvec::SparseMatrix::~SparseMatrix (  void    )  [inline]
 

Definition at line 80 of file sparsematrix.h.

References release(). 

00080 {release();}                      // Destructor

Member Function Documentation

double* matvec::SparseMatrix::a (  void    )  [inline]
 

Definition at line 105 of file sparsematrix.h.

References aav.

Referenced by add_diag(), matvec::GLMM::build_hInv(), close(), matvec::GLMM::contrast(), matvec::Model::get_lms_kp(), getaij(), gibbs_iterate(), initialize_node_list(), insert(), minfilsymelm_node_list(), q(), qTLW(), reset(), matvec::GLMM::residual(), row(), matvec::GLMM::SSQCP(), matvec::Model::uAu_trCA(), matvec::Model::update_mme(), and matvec::Model::vce_gibbs_sampler(). 

00105 {return aav-1;};

void matvec::SparseMatrix::add_diag (  const int    ibeg,
const int    iend,
const double    r
) 
 

Add r to the diagonals from row ibeg to row iend.

Definition at line 731 of file sparsematrix.cpp.

References a(), aav, ia(), iiv, ja(), and jjv.

Referenced by matvec::Model::update_mme(). 

00732 {
00733    int   *ia = iiv-1;
00734    int   *ja = jjv-1;
00735    double *a = aav-1;
00736    int i,j;
00737    for (i=ibeg; i<=iend; i++) {
00738       for (j=ia[i]; j<ia[i+1]; j++) if (ja[j] == i) a[j] += r;
00739    }
00740 }

unsigned matvec::SparseMatrix::close (  void    ) 
 

After close(), no more element can be added.

Definition at line 1060 of file sparsematrix.cpp.

References a(), aav, dim, hsize, matvec::hsort_ija(), ia(), iiv, insertend, ja(), jjv, nonzero, matvec::squeeze(), and matvec::warning().

Referenced by matvec::GLMM::ainv(), matvec::GLMM::save(), matvec::GLMM::setup_mme(), solve(), matvec::Model::uAu_trCA(), matvec::Model::vce_emreml_multi_trait(), and matvec::Model::vce_emreml_single_trait(). 

01061 {                      /* replace iiv by starting elements of rows:*/
01062    if (insertend) return nonzero;
01063    unsigned nrow   = 0;         
01064    unsigned oldrow = 0;
01065    int *ia   = iiv-1;
01066    int *ja   = jjv-1;
01067    double *a = aav-1;
01068    nonzero = squeeze(hsize,ia,ja,a);      /* squeeze out zeros */
01069    hsort_ija(nonzero,ia,ja,a);            /* sort iiv,jjv      */
01070    unsigned k;
01071    for (k=1; k<=nonzero ; k++) {
01072       if (ia[k] != oldrow) {
01073          nrow++;
01074          ia[nrow] = k;
01075          oldrow = ia[k];
01076       }
01077    }
01078    if (nrow != dim) {
01079       throw exception("SparseMatrix::close(): actual size <. expectec size");
01080       // filling diagonals with zeros is one way to avoid this ERROR
01081    }
01082    unsigned i;
01083    ia[dim+1] = nonzero+1;
01084    int zero_pivot = 0;          // zero diagonals ?
01085    for (k=1; k<=dim; k++) {
01086       i = ia[k];
01087       while ( i<ia[k+1] && ja[i] != k ) i++;
01088       if (i < ia[k+1]  && a[i] == 0.0) {
01089          a[i] = 1.0;
01090          zero_pivot++;
01091       }
01092    }
01093    if (zero_pivot) {
01094       warning("zero-diagonals is found in a SparseMatrix, 1.0 has been added");
01095    }
01096    insertend = 1;
01097    //   std::cout << "\n--> A Sparse Matrix has been Set up" << std::endl;
01098    return nonzero;
01099 }

void matvec::SparseMatrix::copyfrom (  const SparseMatrix &    A )  [protected]
 

Get a copy from A

Definition at line 233 of file sparsematrix.cpp.

References aav, dim, elim_done, factor_done, hash_srt, hsize, iiv, info_vec, initial_lij_done, initialize_node_list_done, insertend, jjv, max_nz, nonzero, resize(), solver_name, and srt_hash.

Referenced by operator=(), and SparseMatrix(). 

00234 {
00235    if (this == &A) return;
00236    if (dim != A.dim || max_nz != A.max_nz) resize(A.dim,A.max_nz);
00237    nonzero      = A.nonzero;
00238    hsize        = A.hsize;
00239    info_vec     = A.info_vec;
00240    insertend    = A.insertend;
00241    solver_name  = A.solver_name;
00242    factor_done  = A.factor_done;
00243    elim_done    = A.elim_done;
00244    initial_lij_done = A.initial_lij_done;
00245    initialize_node_list_done = A.initialize_node_list_done;
00246 
00247    memcpy(iiv,A.iiv,sizeof(int)*hsize);
00248    memcpy(jjv,A.jjv,sizeof(int)*hsize);
00249    memcpy(aav,A.aav,sizeof(double)*hsize);
00250    if(hsize){
00251      memcpy(hash_srt,A.hash_srt,sizeof(int)*(hsize+1));
00252      memcpy(srt_hash,A.srt_hash,sizeof(int)*(hsize+1));
00253    }
00254 }

doubleMatrix matvec::SparseMatrix::dense (  void    )  [inline]
 

Definition at line 107 of file sparsematrix.h.

References dim.

Referenced by display(). 

00107 {return dense(1,dim,1,dim);};

doubleMatrix matvec::SparseMatrix::dense (  unsigned    r1,
unsigned    r2,
unsigned    c1,
unsigned    c2
)  const
 

Convert a sub sparse matrix(r1:r2,c1:c2) into Matrix object.

Definition at line 1105 of file sparsematrix.cpp.

Referenced by matvec::GLMM::build_hInv(), and matvec::GLMM::glim(). 

01106 {
01107   doubleMatrix out(r2-r1+1,c2-c1+1);
01108   int   *ia = iiv-1;
01109   int   *ja = jjv-1;
01110   double *a = aav-1;
01111 
01112   unsigned i,j,k;
01113 
01114   for (k=1; k<= hsize; k++) {
01115     if (i = ia[k]) {
01116         j = ja[k];
01117         if (i!=j) {
01118           if (i>=r1 && i<=r2 && j>=c1 && j<=c2 ) {
01119             out(i-r1+1,j-c1+1) = a[k];
01120           }
01121           if (j>=r1 && j<=r2 && i>=c1 && i<=c2 ) {
01122             out(j-r1+1,i-c1+1) = a[k];
01123           }
01124         }
01125         else {
01126           if (i>=r1 && i<=r2) {
01127            out(i-r1+1,i-r1+1) = a[k];
01128           }
01129         }
01130     }
01131   }
01132   return out;
01133 }

void matvec::SparseMatrix::display (  const std::string &    msg = "" )  [inline]
 

Definition at line 127 of file sparsematrix.h.

References dim, and display(). 

00127 {display(1,dim,1,dim,msg);};

void matvec::SparseMatrix::display (  unsigned    r1,
unsigned    r2,
unsigned    c1,
unsigned    c2,
const std::string &    msg = ""
) 
 

Display the content.

Definition at line 1138 of file sparsematrix.cpp.

References dense(), and dim.

Referenced by display(). 

01139 {
01140   if (msg != "") std::cout <<"\n     " << msg << "\n";
01141   if (dim == 0) {
01142     std::cout << "      null sparse matrix     \n\n";
01143     return;
01144   }
01145   std::cout << dense(r1,r2,c1,c2);
01146 }

void matvec::SparseMatrix::display_inv_order (  void    ) 
 

Definition at line 1999 of file sparsematrix.cpp.

References dim, and inv_order. 

01999                                          {
02000   // displays the order of diagonal nodes
02001   unsigned i;
02002   std::cout << "\n" << "inv_order" << "\n";
02003   for (i=1;i<=dim;i++) {
02004     std::cout << inv_order(i) << std::endl;
02005   }
02006   std::cout << std::endl;
02007 }

void matvec::SparseMatrix::display_lij (  void    ) 
 

Definition at line 1975 of file sparsematrix.cpp.

References dim, node_list, and matvec::Vector< NodeStruct >::size(). 

01975                                    {
01976   // displays the lij's
01977   unsigned i,j;
01978   unsigned count;
01979 
01980   for (i=1;i<=dim;i++) {
01981     std::cout <<"\n" <<"node: "<<i<<" ";
01982     count = node_list(i).adj_list.size();
01983     for (j=1;j<=count;j++){
01984       std::cout << node_list(i).a(j) <<" ";
01985     }
01986   }
01987   std::cout << std::endl;
01988 }

void matvec::SparseMatrix::display_node_list (  void    ) 
 

Definition at line 1944 of file sparsematrix.cpp.

References dim, node_list, order, and matvec::Vector< NodeStruct >::size(). 

01944                                          {
01945   // displays the nodes in the adjacent list
01946   // using peeling order numbers
01947   unsigned i,j;
01948   unsigned count;
01949 
01950   for (i=1;i<=dim;i++) {
01951     int orderi = order(i);
01952     std::cout <<"\n" <<"node: " << i <<": ";
01953     
01954     count = node_list(orderi).adj_list.size();
01955     for (j=1;j<=count;j++){
01956       std::cout << node_list(orderi).adj_list(j)<<" ";
01957     }
01958   }
01959   std::cout << std::endl;
01960 
01961 //   unsigned i,j;
01962 //   unsigned count;
01963 //   for (i=1;i<=dim;i++) {
01964 //     std::cout <<"\n" <<"node: "<<i<<" ";
01965 //     count = node_list(i).adj_list.size();
01966 //     for (j=1;j<=count;j++){
01967 //       std::cout << node_list(i).adj_list(j) <<" ";
01968 //     }
01969 //   }
01970 //   std::cout << std::endl;
01971 
01972 
01973 
01974 }

void matvec::SparseMatrix::display_order (  void    ) 
 

Definition at line 1989 of file sparsematrix.cpp.

References dim, and order. 

01989                                      {
01990   // displays the order of diagonal nodes
01991   unsigned i;
01992   std::cout << "\n" << "order" << "\n";
01993   for (i=1;i<=dim;i++) {
01994     std::cout << order(i) << std::endl;
01995   }
01996   std::cout << std::endl;
01997 }

int matvec::SparseMatrix::elim_node (  unsigned    i,
Vector< unsigned >    Aorder
) 
 

Definition at line 2119 of file sparsematrix.cpp.

References matvec::NodeStruct::adj_list, matvec::Vector< unsigned >::begin(), mask, masklist, node_list, matvec::Vector< unsigned >::reserve(), and matvec::Vector< unsigned >::size().

Referenced by minfilsymelm_node_list(). 

02120 {
02121   
02122   // updates elimination graph for eliminating elimnode
02123   
02124   unsigned j,k;
02125   unsigned naborj,nabork,nabor_nabor;
02126   unsigned degnode,degnabor;
02127   unsigned count,*adjpt,*nadjpt,*aadjpt;
02128   unsigned *maskptr,nmask;
02129   //  Vector<unsigned> masklist(dim);
02130   NodeStruct *elimnode, *nabornode;
02131   elimnode = &node_list(Aorder(i));
02132   maskptr=mask.begin();
02133   // loop to update adjacent lists for neigbors of elimnode
02134   degnode = elimnode->adj_list.size();
02135   //  memset(mask.begin(),'\0',dim*sizeof(unsigned)); // SDK
02136 
02137 
02138  
02139   adjpt=elimnode->adj_list.begin();
02140   for(k=1;k<=degnode;k++) {
02141     naborj = *adjpt++;
02142     mask[naborj-1]=1;
02143     masklist[k-1]=naborj-1;
02144   }
02145  
02146   adjpt=elimnode->adj_list.begin();
02147   for (j=1;j<=degnode;j++) { 
02148  
02149     naborj = *adjpt++;
02150     
02151 
02152     nabornode = &node_list(naborj);
02153     nmask=degnode;
02154     degnabor = nabornode->adj_list.size();
02155     nadjpt = nabornode->adj_list.begin();
02156     for (k=0;k<degnabor;k++){
02157       nabor_nabor = *nadjpt++;
02158       if (nabor_nabor != Aorder(i)) { //elim elimnode
02159         if(!mask[nabor_nabor-1]) {  //SDK
02160           masklist[nmask++]=nabor_nabor-1; //SDK
02161         }                      //SDK
02162       }
02163     }
02164     count=nmask-1; //SDK
02165     if (nabornode->adj_list.size() != count) {
02166       nabornode->adj_list.reserve(count);
02167       if (nabornode->adj_list.size() < count) return 0;
02168     } 
02169     int adj_count=0;
02170     nadjpt=nabornode->adj_list.begin();
02171     for(count=0;count<nmask;count++){                  //SDK
02172       if(masklist[count] != (naborj-1)) {
02173         adj_count++;
02174         *nadjpt++ = masklist[count]+1;  //SDK
02175       }
02176     }
02177     if(count != (adj_count+1)) {
02178       std::cout << count << " " << adj_count << "\n" << std::flush;
02179       throw exception("Elim node: not possible");
02180     }//SDK
02181   }
02182   
02183   adjpt=elimnode->adj_list.begin();
02184   for(k=1;k<=degnode;k++) {
02185     naborj = *adjpt++;
02186     mask[naborj-1]=0;
02187   }
02188   return 1;
02189 }

int matvec::SparseMatrix::factor (  void    ) 
 

Definition at line 1667 of file sparsematrix.cpp.

References matvec::NodeStruct::a, matvec::NodeStruct::adj_list, matvec::Vector< NodeStruct >::begin(), matvec::Vector< unsigned >::begin(), matvec::Vector< double >::begin(), matvec::Vector< T >::begin(), diag, dim, factor_done, info_vec, initial_lij(), initial_lij_done, node_list, order, matvec::Vector< T >::reserve(), and matvec::Vector< unsigned >::size().

Referenced by matvec::GLMM::build_hInv(), factorization(), logdet(), matvec::GLMM::restricted_log_likelihood(), and solve(). 

01667                             {
01668   unsigned i,j,nextj,k,count,m;
01669   unsigned coli,colj,rowm,rowk;
01670   unsigned degi,degj;
01671   unsigned startj;
01672   Vector<double> temp;
01673   Vector<unsigned> start, contr_cols;
01674   double lij,diagi;
01675   unsigned *prows;
01676   double   *plij;
01677   NodeStruct *nodes, *node;
01678 
01679   if (!initial_lij_done) {
01680     if(!initial_lij()) return 0;
01681   }
01682   temp.reserve(dim);
01683   start.reserve(dim);
01684   contr_cols.reserve(dim);
01685 
01686 
01687   double   *ptemp        = temp.begin();       ptemp--;
01688   double   *pdiag        = diag.begin();       pdiag--;
01689   unsigned *pstart       = start.begin();      pstart--;
01690   unsigned *porder       = order.begin();      porder--;
01691   unsigned *pcontri_cols = contr_cols.begin(); pcontri_cols--;
01692   double    mult;
01693   unsigned lastpos;
01694   lastpos=0;
01695 
01696   info_vec[0]=0;
01697   info_vec[1]=0;
01698   //for (i=1;i<=dim;i++) {
01699   //  ptemp[i]        = 0.0;
01700   //  pstart[i]       = 0;
01701   //  pcontri_cols[i] = 0;
01702   // }
01703   memset(contr_cols.begin(),'\0',dim*sizeof(unsigned));
01704   memset(start.begin(),'\0',dim*sizeof(unsigned));
01705   memset(temp.begin(),'\0',dim*sizeof(double));
01706   // loop for calculating column i of lij
01707   nodes = node_list.begin(); nodes--;
01708   for (i=1;i<=dim;i++){
01709     coli = porder[i];
01710     j = pcontri_cols[i];
01711     count = 0;
01712     // loop for getting contributions from columns contributing to column i
01713     // contri_cols is used to point to these columns
01714     while (j) {
01715       startj = pstart[j];
01716       colj  = porder[j];
01717       node  = &nodes[colj];
01718       degj  = node->adj_list.size();
01719       nextj = pcontri_cols[j];
01720       prows = node->adj_list.begin();
01721       plij  = node->a.begin();
01722       lij = plij[startj];
01723       // update pstart[j] for next time
01724       pstart[j] = startj+1;
01725       ptemp[prows[startj++]]+=lij*lij;
01726       // accumulate contributions from col j to col i in temp
01727       for (k=startj;k <degj ;k++) {
01728         rowk=prows[k];
01729         ptemp[rowk]  =ptemp[rowk]+ lij*plij[k] ;
01730       }
01731     
01732       
01733       k = startj;
01734       if (k<degj) {
01735         rowk = prows[k];
01736         pcontri_cols[j] = pcontri_cols[rowk];
01737         pcontri_cols[rowk] = j;
01738        }
01739       j = nextj;
01740     }
01741     
01742     node  = &nodes[coli];
01743     prows = node->adj_list.begin();
01744     plij  = node->a.begin();
01745     diagi = pdiag[coli] - ptemp[i];
01746     degi = node->adj_list.size();
01747     //         std::cout << "No of contributions to row: " << count << std::endl;
01748     //         std::cout << "No of adj elements       : " << degi   << std::endl;
01749     //         std::cout << "diagonal element: " << coli << " with value: " << pdiag[coli] << std::endl;
01750     //         std::cout << "row             : " << i    << " with value: " << ptemp[i] << std::endl;
01751     //         std::cout << "check psd:" << diagi << std::endl;
01752     //         std::cout << std::endl;
01753     if (diagi > std::max(1.e-8,1.e-10*pdiag[coli])) { 
01754       info_vec[0]+=std::log(diagi);
01755       info_vec[1]++;
01756       //std::cout << " rank " << info_vec[1] ;
01757       diagi = std::sqrt(diagi);
01758       pdiag[coli] = diagi;
01759       mult = 1/diagi;
01760       // use contributions in temp to update col i
01761       for (k=0;k<degi;k++){
01762         rowk = prows[k];
01763         plij[k]-=ptemp[rowk];
01764         plij[k]*=mult;
01765         ptemp[rowk] = 0.0;
01766       }
01767       
01768       if (degi) {
01769         rowk = prows[0];
01770         pcontri_cols[i] = pcontri_cols[rowk];
01771         pcontri_cols[rowk] = i;
01772       }
01773       
01774     }
01775     else if (diagi > -std::max(1.e-10,1.e-10*pdiag[coli])) {
01776       pdiag[coli] = 0;
01777       degi = node->adj_list.size();
01778       for (k=0;k<degi;k++){
01779         rowk = prows[k];
01780         plij[k] = 0.0;
01781         ptemp[rowk] = 0.0;
01782       }
01783     }
01784     else if (diagi < 0) {
01785       throw exception("SparseMatrix::factor matrix is not positive semi-definite");
01786     }
01787   }
01788 
01789   factor_done = 1;
01790   return 1;
01791 }

unsigned matvec::SparseMatrix::factorization (  const int    path = 5,
const double    tol = 1.0e-10
) 
 

Factorization.

Parameters:
path  an integer path specification; values and their meanings are -
  • 4 perform SSF
  • 5 perform SSF and SNF
  • 6 perform SNF only (isp/rsp is assumed to have been set up in an earlier call to sdrv (for SSF))

Definition at line 868 of file sparsematrix.cpp.

References factor(), initial_lij(), initialize_node_list(), irank(), and minfilsymelm_node_list().

Referenced by matvec::Model::graph_log_likelihood_peeling(), matvec::Model::log_likelihood_peeling(), and matvec::Model::vce_emreml_multi_trait(). 

00869 {
00870    if (path < 4) throw exception("SparseMatrix::factorization(path): path >= 4 is required");
00871 
00872    if (path==4) {
00873      initialize_node_list();
00874      minfilsymelm_node_list();
00875      initial_lij();
00876      factor();
00877    }
00878    if (path==6){
00879      factor();
00880    }
00881    return irank();
00882 }

double matvec::SparseMatrix::getaij (  const int    ii,
const int    jj
) 
 

Return the element at (i,j) position.

Definition at line 977 of file sparsematrix.cpp.

References a(), aav, dim, hsize, ia(), iiv, ja(), and jjv.

Referenced by initial_lij(), operator()(), matvec::GLMM::residual(), and matvec::GLMM::SSQCP(). 

00978 {
00979   int i,j;
00980   i=ii;
00981   j=jj;
00982   
00983    if (dim == 0) throw exception("SparseMatrix is empty");
00984    if (i>j) {
00985      i=jj;
00986      j=ii;
00987    }
00988 //throw exception("SparseMatrix::insert(i,j,a): i>j, only upper triangular allowed");
00989    if (j>dim) {
00990      //std::cout << "getaij "<< i<<","<<j;
00991     throw exception("SparseMatrix::getaij(i,j): range error"); 
00992   }
00993 
00994    int *ia   = iiv-1;
00995    int *ja   = jjv-1;
00996    double *a = aav-1;
00997    unsigned new_ill, ill;
00998    ill =  1 + (433*i+53*j) % hsize;
00999    for (unsigned k=1; k<=200; k++) {
01000       if (ia[ill] == 0) {       // aij = 0
01001         return(0.0);
01002       }
01003       else if (ia[ill] == i && ja[ill] == j) {
01004          return(a[ill]);
01005       }
01006       else {
01007          new_ill = (ill+640) % hsize + 1;
01008          ill = (new_ill == ill) ? (ill+1)% hsize + 1 : new_ill;
01009       }
01010    }
01011    throw exception("SparseMatrix::getaij This should not happen! Probably a bug :-( ");
01012    return(0.0);
01013 }

void matvec::SparseMatrix::gibbs_iterate (  double *    x,
const double *    r,
double *    w,
const double    se = 1.0
) 
 

Definition at line 1175 of file sparsematrix.cpp.

References a(), aav, dim, ia(), iiv, ja(), and jjv.

Referenced by matvec::Model::genotype_dist_gibbs(), matvec::Model::genotypic_value_cat(), and matvec::Model::genotypic_value_gibbs(). 

01177 {
01178    /////////////////////////////////////////////////////////////////
01179    // REQUIREMENTS: x,r, and w must have enough spaces
01180    // se = sqrt(residual variance), if MMe is built with variance
01181    // ratio, then se is necessary, otherwise set 1 which is default
01182    //////////////////////////////////////////////////////////////////
01183    double *sample = x-1;
01184    double *adj_rhs = w-1;
01185    double *a = aav-1;
01186    int *ia = iiv-1;
01187    int *ja = jjv-1;
01188    int i,j,jj,k,ke;
01189    double mean,diagon=0.0;
01190 
01191    memcpy(w,r,sizeof(double)*dim);
01192    for (i=1; i<=dim; i++) {
01193       mean = adj_rhs[i];
01194       k = ia[i];  ke = ia[i+1];
01195       for (j=k; j<ke; j++) {
01196          jj = ja[j];
01197          if (jj == i) {diagon = a[j];}  // diagonal may not be the first elements
01198          else         {mean -= a[j]*sample[jj];}
01199       }
01200       mean /= diagon;
01201 //      mean += se*snorm()/std::sqrt(diagon);   // variance = 1/diagon
01202       sample[i] = mean;                // adjust the rhs, use all coef of row i
01203       for (j=k; j<ke; j++)  adj_rhs[ja[j]] -= a[j]*mean;
01204    }
01205 }

int* matvec::SparseMatrix::ia (  void    )  [inline]
 

Definition at line 117 of file sparsematrix.h.

References iiv.

Referenced by add_diag(), matvec::GLMM::build_hInv(), close(), matvec::GLMM::contrast(), matvec::Model::get_lms_kp(), getaij(), gibbs_iterate(), initialize_node_list(), insert(), q(), qTLW(), reset(), matvec::GLMM::residual(), row(), matvec::GLMM::SSQCP(), matvec::Model::uAu_trCA(), matvec::Model::update_mme(), and matvec::Model::vce_gibbs_sampler(). 

00117 {return iiv-1;};

int matvec::SparseMatrix::initial_lij (  void    ) 
 

Definition at line 1626 of file sparsematrix.cpp.

References matvec::NodeStruct::a, matvec::NodeStruct::adj_list, matvec::Vector< double >::begin(), matvec::Vector< unsigned >::begin(), matvec::Vector< NodeStruct >::begin(), diag, dim, elim_done, getaij(), initial_lij_done, minfilsymelm_node_list(), node_list, order, row(), and matvec::Vector< unsigned >::size().

Referenced by factor(), factorization(), and reorder(). 

01627 {
01628   if (insertend) throw exception("SparseMatrix::initial_lij: cannot initial_lij after closing");
01629   unsigned i,j;
01630   unsigned row,col;
01631   unsigned deg,nabor;
01632   double aij;
01633   double *paij;
01634   unsigned *nabors;
01635   NodeStruct *nodes, *node;
01636 
01637   if (!elim_done) {
01638     if (!minfilsymelm_node_list()) return 0;
01639   }
01640 
01641   nodes = node_list.begin();
01642   for (i=0;i<dim;i++) {
01643     row = i+1;
01644     diag[i] = getaij(row,row);
01645     node = &nodes[i];
01646     nabors = node->adj_list.begin();
01647     deg    = node->adj_list.size();
01648 
01649     for (j=0;j<deg;j++){
01650       nabor = nabors[j];
01651       col = order(nabor);
01652       if (col < row) {
01653         aij = getaij(col,row);
01654       }
01655       else{
01656         aij = getaij(row,col);
01657       }
01658       paij = node->a.begin();
01659       paij[j] = aij;
01660     }
01661   }
01662   initial_lij_done = 1;
01663   return 1;
01664 }

int matvec::SparseMatrix::initialize_node_list (  void    ) 
 

Definition at line 1298 of file sparsematrix.cpp.

References a(), aav, dim, ia(), iiv, initialize_node_list_done, ja(), jjv, mask, node_list, nonzero, matvec::Vector< NodeStruct >::reserve(), matvec::Vector< unsigned >::resize(), row(), and srt_hash.

Referenced by factorization(), minfilsymelm_node_list(), and reorder(). 

01298                                            {
01299   int *ia   = iiv-1;
01300   int *ja   = jjv-1;
01301   double *a = aav-1;
01302   unsigned i,j;
01303   unsigned n_rows,row,col;
01304   unsigned count;
01305 
01306   // I am going to count number of nabors for each node
01307   // I will use "mask" to do the counting!
01308 
01309   mask.resize(dim);
01310 
01311   node_list.reserve(dim);
01312   // do the counting
01313   //SDK
01314   int isrt;
01315   for(isrt=0;isrt<nonzero;isrt++) {
01316     i=srt_hash[isrt];
01317     //  for (i=1;i<=hsize;i++){
01318     //SDK
01319     if (!ia[i]) continue; //SDK
01320     row = ia[i];
01321     col = ja[i];
01322     if(row != col){
01323       mask(row) = mask(row) + 1;
01324       mask(col) = mask(col) + 1;
01325     }
01326   }
01327     // now resize the adjacent list for each node
01328   for (i=1;i<=dim;i++){
01329     node_list(i).adj_list.reserve(mask(i));
01330   }
01331   // now I am going to use mask to mark the current position in the adj_list for each node
01332   // start at position 1
01333   for (i=1;i<=dim;i++){
01334     mask(i) = 1;
01335   }
01336     // now do the filling
01337   //SDK
01338   for(isrt=0;isrt<nonzero;isrt++) {
01339     i=srt_hash[isrt];
01340     //  for (i=1;i<=hsize;i++){
01341     //SDK
01342     if (!ia[i]) continue; //SDK
01343     row = ia[i];
01344     col = ja[i];
01345     if(row!= col){
01346       node_list(row).adj_list(mask(row)) = col;
01347       node_list(col).adj_list(mask(col)) = row;
01348       mask(row) = mask(row) + 1;
01349       mask(col) = mask(col) + 1;
01350     }
01351   }
01352   initialize_node_list_done = 1;
01353   return 1;
01354 }

int matvec::SparseMatrix::insert (  const int    i,
const int    j,
const double    value
) 
 

Add value to element (i,j).

Note that Indexing is from 1 instead of 0

Definition at line 927 of file sparsematrix.cpp.

References a(), aav, dim, matvec::Session::epsilon, hash_srt, hsize, ia(), iiv, ja(), jjv, nonzero, matvec::SESSION, and srt_hash.

Referenced by matvec::Model::add_Ag(), matvec::GLMM::add_Ag(), matvec::GLMM::add_Ag_old(), matvec::GLMM::add_AgSand(), matvec::Model::add_G_1_single_trait(), matvec::Population::add_Ga_inv(), matvec::Model::add_Ig(), matvec::GLMM::add_IgSand(), matvec::GLMM::ainv(), matvec::GLMM::build_hInv(), resize(), matvec::Population::setup_m_ww(), matvec::Model::setup_ww(), matvec::GLMM::setup_ww(), matvec::Model::setup_ww_single_trait(), SparseMatrix(), and matvec::Model::vce_emreml_multi_trait(). 

00928 {
00929   int ii,jj; 
00930   ii = i;    
00931   jj = j;    
00932   
00933   if (dim == 0) throw exception("SparseMatrix::insert(): SparseMatrix is empty");
00934   if (i<=0 || j<=0 || j>dim) {
00935     //std::cout << "insert "<< i<<","<<j;
00936     throw exception("SparseMatrix::insert(i,j,a): range error"); 
00937   }
00938   if (i>j) {
00939     // symmetry is assumed, and only upper triangular elements saved
00940     ii = j;
00941     jj = i;
00942   }
00943   //    if (i>j) throw exception("SparseMatrix::insert(i,j,a): i>j, only upper triangular allowed");
00944 
00945   if (fabs(value) <= SESSION.epsilon && ii != jj) return(1);
00946   int *ia   = iiv-1;
00947   int *ja   = jjv-1;
00948   double *a = aav-1;
00949   unsigned new_ill, ill;
00950   ill =  1 + (433*ii+53*jj) % hsize;
00951   for (unsigned k=1; k<=200; k++) {
00952     if (ia[ill] == 0) { // first time store
00953       srt_hash[nonzero]=ill;
00954       hash_srt[ill]=nonzero;
00955       nonzero++;  //SDK
00956       ia[ill]=ii;
00957       ja[ill]=jj;
00958       a[ill]=value;
00959       return(1);
00960     }
00961     else if (ia[ill] == ii && ja[ill] == jj) {
00962       a[ill] += value;
00963       return(1);
00964     }
00965     else {
00966       new_ill = (ill+640) % hsize + 1;
00967       ill = (new_ill == ill) ? (ill+1)% hsize + 1 : new_ill;
00968     }
00969   }
00970   throw exception("SparseMatrix::insert(), max_nz is too small");
00971   return(0);
00972 }

unsigned matvec::SparseMatrix::irank (  void    ) 
 

Return the rank of the obhect.

Definition at line 897 of file sparsematrix.cpp.

References logdet(), logdet_done, and rank.

Referenced by factorization(), logdet(), and matvec::Model::vce_emreml_single_trait(). 

00897                              {
00898   if (!logdet_done) {
00899     logdet(rank);
00900   }
00901   return rank;
00902 }

int* matvec::SparseMatrix::ja (  void    )  [inline]
 

Definition at line 118 of file sparsematrix.h.

References jjv.

Referenced by add_diag(), matvec::GLMM::build_hInv(), close(), matvec::GLMM::contrast(), matvec::Model::get_lms_kp(), getaij(), gibbs_iterate(), initialize_node_list(), insert(), q(), qTLW(), reset(), matvec::GLMM::residual(), row(), matvec::GLMM::SSQCP(), matvec::Model::uAu_trCA(), matvec::Model::update_mme(), and matvec::Model::vce_gibbs_sampler(). 

00118 {return jjv-1;};

double matvec::SparseMatrix::logdet (  void    ) 
 

Returns log of its determinant.

Definition at line 887 of file sparsematrix.cpp.

References logdet_done, logdeterm, and rank.

Referenced by irank(), and logdet(). 

00887                             {
00888   if (!logdet_done) {
00889     logdet(rank);
00890   }
00891   return logdeterm;
00892 }

double matvec::SparseMatrix::logdet (  unsigned *    irank )  [inline]
 

Definition at line 110 of file sparsematrix.h.

References irank(), and logdet(). 

00110 {return logdet(*irank);};

double matvec::SparseMatrix::logdet (  unsigned &    irank ) 
 

Definition at line 904 of file sparsematrix.cpp.

References diag, dim, factor(), factor_done, irank(), and logdeterm.

Referenced by matvec::GLMM::ainv(), matvec::Model::anova(), matvec::GLMM::log_likelihood(), and matvec::Model::vce_emreml_single_trait(). 

00904                                             {
00905    if (dim == 0) throw exception("SparseMatrix::logdet(): SparseMatrix is empty");
00906   unsigned i;
00907   double res=0.0;
00908    if (!factor_done) {
00909      if(!factor()) return 0.0;
00910    }
00911    irank = 0;
00912    for (i=1;i<=dim;i++){
00913      if (diag(i) > 0.0000000001) {
00914        irank++;
00915        res += std::log(diag(i));
00916      }
00917    }
00918    logdeterm = 2*res;
00919    return (logdeterm);
00920 }

int matvec::SparseMatrix::minfilsymelm_node_list (  void    ) 
 

Definition at line 1356 of file sparsematrix.cpp.

References a(), adjncy, matvec::Vector< NodeStruct >::begin(), matvec::cmp_func, matvec::comp(), diag, dim, elim_done, elim_node(), matvec::idxtype, initialize_node_list(), initialize_node_list_done, inv_order, mask, masklist, node_list, nonzero, order, matvec::Vector< double >::reserve(), matvec::Vector< NodeStruct >::reserve(), matvec::Vector< unsigned >::reserve(), matvec::Vector< unsigned >::resize(), matvec::Vector< T >::resize(), matvec::Vector< NodeStruct >::size(), and xadj.

Referenced by factorization(), initial_lij(), and reorder(). 

01357 {
01358   unsigned i,j;
01359   unsigned orderj,inor;
01360   unsigned mindeg=dim,deg,oldmindeg=0;
01361   Vector<unsigned> same;
01362   unsigned nsame;
01363   unsigned minpos,elimpos,temp;
01364   unsigned nabor;
01365   unsigned count;
01366   unsigned *num;
01367   //std::cout << "entering minfilsymelm_node_list" << std::endl;
01368 
01369 
01370   if (!initialize_node_list_done) {
01371     if (!initialize_node_list()) return 0;
01372   }
01373 
01374   order.reserve(dim);
01375   inv_order.reserve(dim);
01376   same.resize(dim);
01377   // initialize order
01378   for (i=1;i<=dim;i++) {
01379     order(i) = i;
01380   }
01381 
01382 
01383   if(xadj) delete [] xadj;
01384   if(adjncy) delete [] adjncy;
01385   xadj=new idxtype [dim+1];
01386   adjncy=new idxtype [2*(nonzero-dim)];
01387   xadj[0]=1;
01388   idxtype *xoff,*aoff;
01389   xoff=xadj;xoff--;
01390   aoff=adjncy;aoff--;
01391   idxtype idcnt=1;
01392   for(i=1;i<=dim;i++) {
01393       deg = node_list(i).adj_list.size();
01394       xoff[i+1]=xoff[i]+deg; 
01395       for(j=1;j<=deg;j++,idcnt++) {
01396         aoff[idcnt]= node_list(i).adj_list(j);
01397       } 
01398   }
01399   masklist.reserve(dim);
01400   mask.resize(dim);
01401   
01402 #ifdef DO_MINORD
01403   // loop for elimination
01404   Vector<unsigned> degpt,degnxt,deglst;
01405   degpt.resize(dim+1);
01406   degnxt.resize(dim+1);
01407   deglst.resize(dim+1);
01408   unsigned tmp1,tmp2,tmp3,k;
01409   for(i=1;i<=dim;i++){
01410     deg=node_list(i).adj_list.size();
01411     tmp1=degpt[deg];
01412     degpt[deg]=i;
01413     degnxt[i]=tmp1;
01414     if(tmp1) deglst[tmp1]=i;
01415   }
01416   //  cout << "degpt\n" << degpt; 
01417   //  cout << "degnxt\n" << degnxt; 
01418   //  cout << "deglst\n" << deglst; 
01419   mindeg=1;
01420   for (i=1;i<=dim;i++) {
01421     // loop to look through all remaining nodes to update mindeg
01422     //mindeg = dim+1;
01423     nsame=0;
01424     if(mindeg) mindeg--;
01425     for(k=mindeg;degpt[k]==0&&k<=dim;k++);
01426     minpos=degpt[k];
01427     //if(k<mindeg) cout << "Mindeg " << mindeg <<" " << k << std::endl;
01428     mindeg=k;
01429     /*
01430     for (j=i;j<=dim;j++) {
01431       orderj = order(j);
01432       deg = node_list(orderj).adj_list.size();
01433       if (deg < mindeg) {
01434         mindeg = deg;
01435         minpos = j;
01436       }
01437     }
01438     */
01439     
01440     //    cout << "Mindeg " << mindeg << std::endl;
01441     // update order vector
01442     for(j=0;j<mindeg;j++) {
01443       //orderj = order(i);
01444       k=node_list(minpos).adj_list[j];
01445       same[j]=node_list(k).adj_list.size();
01446     }
01447     temp = order(i);
01448     order(i)=minpos;
01449     //order(i) = order(minpos);
01450     //order(minpos) = temp;
01451     // now eliminate node at minpos
01452     //    std::cout << i <<" eliminating row at " << minpos<< " " << std::endl;
01453     if (!elim_node(i,order)) return 0;
01454 
01455     {
01456       k=minpos;
01457       unsigned a,b,c;
01458       a=deglst[k];
01459       c=degnxt[k];
01460       deglst[c]=0;
01461       degpt[mindeg]=c;
01462       if(a !=0) std::cout << "Err " << a << " " << k << std::endl;
01463 
01464     }
01465     for(j=0;j<mindeg;j++) {
01466       k=node_list(minpos).adj_list[j];
01467       deg=node_list(k).adj_list.size();
01468       if(same[j]!=deg){
01469         unsigned a,b,c;
01470         b=degpt[deg];
01471         c=degnxt[k];
01472         a=deglst[k];
01473         //if(a == minpos)a=0;
01474         degpt[deg]=k;
01475         degnxt[k]=b;
01476         deglst[k]=0;
01477         if(b) deglst[b]=k;
01478         if(a) degnxt[a]=c;
01479         if(c) deglst[c]=a;
01480         if(a==0) degpt[same[j]]=c;
01481 
01482       }
01483     }
01484     //  cout << "degpt\n" << degpt; 
01485     //  cout << "degnxt\n" << degnxt; 
01486     //  cout << "deglst\n" << deglst; 
01487   }
01488 #endif
01489   for (i = 1;i<=dim;i++) {
01490     // loop to get inv_order from order
01491 #ifndef DO_MINORD
01492     if (!elim_node(i,order)) return 0;
01493 #endif
01494     inor = order(i);
01495     inv_order(inor) = i;
01496   }
01497   
01498   //display_order();
01499   //display_inv_order();
01500   // renumber and rearrange adj_lists according to order of elimination
01501   for (i=1;i<=dim;i++) {
01502     deg = node_list(i).adj_list.size();
01503     // put nabors into mask according to elimination order
01504     for (j=1;j<=deg;j++){
01505       nabor = node_list(i).adj_list(j);
01506       elimpos = inv_order(nabor);
01507       node_list(i).adj_list(j) = elimpos;
01508     }
01509     num = node_list(i).adj_list.begin();
01510     qsort(num,deg,sizeof(unsigned), (cmp_func) comp);
01511   }
01512   // get memory for lij and diagonals
01513   unsigned totsize=0;
01514   for (i=1;i<=dim;i++) {
01515     deg = node_list(i).adj_list.size(); 
01516     node_list(i).a.reserve(deg);
01517     totsize+=deg;
01518     //    std::cout << " i " << i << ":";
01519     //    num = node_list(i).adj_list.begin();
01520     //    for(j=0;j<deg;j++) std::cout << " " << num[j];
01521     //    std::cout << "\n";
01522     if (node_list(i).a.size() < 0) return 0;
01523   }
01524   //  std::cout << "Total : " << totsize << "\n";
01525   diag.reserve(dim);
01526   elim_done = 1;
01527   //std::cout << "exiting minfilsymelm_node_list" << std::endl;
01528   return 1;
01529 }

void matvec::SparseMatrix::mv (  const double *    v,
const unsigned    n,
double *    result
) 
 

Definition at line 436 of file sparsematrix.cpp.

References aav, dim, iiv, jjv, nonzero, and srt_hash. 

00436                                                                        {
00437    /////////////////////////////////////////////////////////////////////
00438    // REQUIREMENTS:  both v and result must have n(dim) elements (space)
00439    // modified 8/4/1998
00440    //////////////////////////////////////////////////////////////////////
00441 
00442   if (dim != n) throw exception("SparseMatrix::mv(): not conformable");
00443   memset(result,'\0',sizeof(double)*n);
00444   int *A_ia   = iiv-1;   // offset by 1
00445   int *A_ja   = jjv-1;
00446   double *A_a = aav-1;
00447   const double *v_ve = v-1;
00448   unsigned i,j,k;
00449   double x;
00450 
00451   //  result.resize(n,0.0);
00452 
00453   result--;      // offset by 1
00454   //SDK
00455   int ksrt;
00456   for(ksrt=0;ksrt<nonzero;ksrt++) {
00457     k=srt_hash[ksrt];
00458     // for (k=1; k<= hsize; k++) {
00459     //SDK
00460     if (i = A_ia[k]) {
00461         j = A_ja[k];
00462         x = A_a[k];
00463         result[i] = result[i] + x*v_ve[j];
00464         if (j>i) {
00465           result[j] = result[j] + x*v_ve[i];
00466         }
00467     }
00468   }
00469   result++; // offset back
00470 }

void matvec::SparseMatrix::mv (  const Vector< double > &    v,
Vector< double > &    result
) 
 

Definition at line 428 of file sparsematrix.cpp.

References matvec::Vector< T >::begin(), dim, matvec::Vector< T >::reserve(), and matvec::Vector< T >::size().

Referenced by matvec::GLMM::contrast(). 

00429 {
00430    if (dim != v.size()) throw exception("SparseMatrix::mv(): sizenot conformable");
00431    if (v.begin() == result.begin()) throw exception("SparseMatrix::mv(): can't be done in situ");
00432    result.reserve(dim);
00433    this->mv(v.begin(),dim,result.begin());
00434 }

unsigned matvec::SparseMatrix::num_rows (  void    )  const [inline]
 

Definition at line 101 of file sparsematrix.h.

References dim.

Referenced by matvec::GLMM::contrast(), matvec::GLMM::residual(), and matvec::GLMM::SSQCP(). 

00101 {return dim;};

unsigned matvec::SparseMatrix::nz (  void    )  const [inline]
 

Definition at line 102 of file sparsematrix.h.

References nonzero.

Referenced by matvec::Model::blup(), matvec::Model::get_blupsol(), matvec::Model::info(), matvec::GLMM::info(), matvec::Model::setup_mme(), and matvec::GLMM::setup_mme(). 

00102 {return nonzero;};

double matvec::SparseMatrix::operator() (  const int    i,
const int    j
) 
 

Retrieve element (i,j).

Definition at line 259 of file sparsematrix.cpp.

References dim, and getaij(). 

00259                                                        {
00260 
00261   if (i>=1 && i<=dim && j>=1 && j<=dim) {
00262     return getaij(i,j);
00263   }
00264   else {
00265     throw exception("SparseMatrix(): subscript out of range");
00266   }
00267   return 0.0;
00268 }

const SparseMatrix & matvec::SparseMatrix::operator= (  const SparseMatrix &    A ) 
 

Assignment operator.

Definition at line 273 of file sparsematrix.cpp.

References copyfrom(). 

00274 {
00275   copyfrom(A);
00276   return *this;
00277 }

double matvec::SparseMatrix::q (  const double *    x,
const double *    y,
const int    i1,
const int    i2
) 
 

Calculates the quadratic form x*A*y from i1'th row to i2'th row.

SparseMatrix A is half storage (upper or lower or mixed) REQUIREMENTS: x and y must have at least this->dim elements x[0] and y[0] are the first elements.

Definition at line 591 of file sparsematrix.cpp.

References a(), aav, dim, ia(), iiv, ja(), jjv, nonzero, and srt_hash. 

00592 {
00593   if (dim == 0) throw exception("SparseMatrix::q(): SparseMatrix is empty");
00594   if (i1<=0 || i1>dim || i2<=0 || i2>dim || i1>i2) throw exception("SparseMatrix::q(): out of range");
00595 
00596   int   *ia = iiv-1;
00597   int   *ja = jjv-1;
00598   double *a = aav-1;
00599   double tmp,qaq ;
00600   const double *xx,*yy;
00601   xx = x-1;  yy = y-1;
00602   unsigned i,j,k;
00603   //SDK
00604   int ksrt;
00605   for (qaq = 0.0, ksrt=0; ksrt< nonzero; ksrt++) {
00606     k=srt_hash[ksrt];
00607     //  for (qaq = 0.0, k=1; k<= hsize; k++) {
00608     //SDK
00609     if (i = ia[k]) {
00610         j = ja[k];
00611         if (i == j && i >= i1 && i <= i2 ) {
00612           tmp = xx[i]*a[k]*yy[i];
00613           qaq += tmp;
00614         }
00615         else if (i >= i1 && i <= i2 && j >= i1 && j <= i2) {
00616           tmp = xx[i]*a[k]*yy[j] + xx[j]*a[k]*yy[i];
00617           qaq += tmp;
00618         }
00619     }
00620   }
00621   return qaq;
00622 }

double matvec::SparseMatrix::q (  const double *    x,
const double *    y
) 
 

Calculates the quadratic form x*A*y.

SparseMatrix A is half storage (upper or lower or mixed) REQUIREMENTS: x and y must have at least this->dim elements x[0] and y[0] are the first elements

Definition at line 552 of file sparsematrix.cpp.

References a(), aav, dim, ia(), iiv, ja(), jjv, nonzero, and srt_hash. 

00553 {
00554   if (dim == 0) exception("SparseMatrix is empty");
00555 
00556   int   *ia = iiv-1;
00557   int   *ja = jjv-1;
00558   double *a = aav-1;
00559   double tmp,qaq ;
00560   const double *xx,*yy;
00561   xx = x-1;  yy = y-1;
00562   unsigned i,j,k;
00563   //SDK
00564   int ksrt;
00565   for(qaq=0.0,ksrt=0;ksrt<nonzero;ksrt++) {
00566     k=srt_hash[ksrt];
00567     //  for (qaq = 0.0, k=1; k<=hsize; k++) {
00568     //SDK
00569     if (i = ia[k]) {
00570         j = ja[k];
00571         if (i == j) {
00572           tmp = xx[i]*a[k]*yy[i];
00573           qaq += tmp;
00574         }
00575         else {
00576           tmp = xx[i]*a[k]*yy[j] + xx[j]*a[k]*yy[i];
00577           qaq += tmp;
00578         }
00579     }
00580   }
00581   return qaq;
00582 }

double matvec::SparseMatrix::q (  const Vector< double > &    v1,
const Vector< double > &    v2
) 
 

Calculates the product v1*A*v2.

SparseMatrix A is half storage only(upper,lower or mixed)

Definition at line 538 of file sparsematrix.cpp.

References matvec::Vector< T >::begin(), dim, and matvec::Vector< T >::size().

Referenced by matvec::Model::graph_log_likelihood_peeling(), matvec::Model::log_likelihood_peeling(), qTLW(), matvec::Model::vce_emreml_multi_trait(), and matvec::Model::vce_gibbs_sampler(). 

00539 {
00540    if (dim == 0) throw exception("SparseMatrix is empty");
00541    if (dim != v1.size() || dim != v2.size()) throw exception("SparseMatrix::q(v1,v2): size not conformable");
00542    return q(v1.begin(),v2.begin());
00543 }

double matvec::SparseMatrix::qTLW (  const double *    x,
const double *    y,
const int    i1,
const int    i2
) 
 

Calculates the quadratic form x*A*y from i1'th row to i2'th row.

SparseMatrix A is half storage (upper or lower or mixed) REQUIREMENTS: x and y must have at least this->dim elements x[0] and y[0] are the first elements

Definition at line 686 of file sparsematrix.cpp.

References a(), aav, dim, ia(), iiv, insertend, ja(), and jjv. 

00688 {
00689    if (dim == 0) throw exception("SparseMatrix::qTLW(): SparseMatrix is empty");
00690    if (!insertend) throw exception("SparseMatrix::qTLW(): call SparseMatrix::close() first");
00691    if (i1<=0 || i1>dim || i2<=0 || i2>dim || i1>i2) throw exception("SparseMatrix::qTLW(), out of range");
00692 
00693    int   *ia = iiv-1;
00694    int   *ja = jjv-1;
00695    double qaq,tmp, *a = aav-1;
00696    const double *xx,*yy;
00697    xx = x-1;  yy = y-1;
00698    int i,j,jj,je;
00699    if (x == y) {
00700       for (qaq=0.0,i=i1; i<=i2; i++) {
00701          je = ia[i+1];
00702          for (j=ia[i]; j<je; j++) {
00703             jj = ja[j];
00704             if (jj >= i1 && jj <= i2) {     // this is necessary since M may not
00705                tmp = xx[i]*a[j]*yy[jj];     // be strictly upper or lower
00706                qaq += tmp;
00707                if (jj != i) qaq += tmp;
00708             }
00709          }
00710       }
00711    }
00712    else {
00713       for (qaq=0.0,i=i1; i<=i2; i++) {
00714          je = ia[i+1];
00715          for (j=ia[i]; j<je; j++) {
00716             jj = ja[j];
00717             if (jj >= i1 && jj <= i2) {     // this is necessary since M may not
00718                tmp = xx[i]*a[j]*yy[jj];     // be strictly upper or lower
00719                qaq += tmp;
00720                if (jj != i) qaq += xx[jj]*a[j]*yy[i];
00721             }
00722          }
00723       }
00724    }
00725    return qaq;
00726 }

double matvec::SparseMatrix::qTLW (  const double *    x,
const double *    y
) 
 

Calculates the quadratic form x*A*y

SparseMatrix A is half storage (upper or lower or mixed) REQUIREMENTS: x and y must have at least this->dim elements x[0] and y[0] are the first elements

Definition at line 643 of file sparsematrix.cpp.

References a(), aav, dim, ia(), iiv, insertend, ja(), and jjv. 

00644 {
00645    if (dim ==0) throw exception("SparseMatrix::qTLW(): SparseMatrix is empty");
00646    if (!insertend) throw exception("SparseMatrix::qTLW(): first call SparseMatrix::close() first");
00647 
00648    int   *ia = iiv-1;
00649    int   *ja = jjv-1;
00650    double qaq,tmp, *a = aav-1;
00651    const double *xx,*yy;
00652    xx = x-1;  yy = y-1;
00653    int i,j,je,jj;
00654    if (x == y) {
00655       for (qaq=0.0,i=1; i<=dim; i++) {
00656          je = ia[i+1];
00657          for (j=ia[i]; j<je; j++) {
00658             jj = ja[j];
00659             tmp = xx[i]*a[j]*yy[jj];
00660             qaq += tmp;
00661             if (jj != i) qaq += tmp;
00662          }
00663       }
00664    }
00665    else {
00666       for (qaq=0.0,i=1; i<=dim; i++) {
00667          je = ia[i+1];
00668          for (j=ia[i]; j<je; j++) {
00669             jj = ja[j];
00670             tmp = xx[i]*a[j]*yy[jj];
00671             qaq += tmp;
00672             if (jj != i) qaq +=  xx[jj]*a[j]*yy[i];
00673          }
00674       }
00675    }
00676    return qaq;
00677 }

double matvec::SparseMatrix::qTLW (  const Vector< double > &    v1,
const Vector< double > &    v2
) 
 

Calculates the quadratic form v1*A*v2.

SparseMatrix A is half storage only(upper,lower or mixed)

Definition at line 629 of file sparsematrix.cpp.

References matvec::Vector< T >::begin(), dim, q(), and matvec::Vector< T >::size().

Referenced by matvec::Model::uAu_trCA(). 

00630 {
00631    if (dim == 0) throw exception("SparseMatrix::qTLW(): it is empty");
00632    if (dim != v1.size() || dim != v2.size()) throw exception("SparseMatrix::qTLW(): size not conformable");
00633    return q(v1.begin(),v2.begin());
00634 }

void matvec::SparseMatrix::release (  void    ) 
 

Definition at line 503 of file sparsematrix.cpp.

References aav, adjncy, dim, hash_srt, hsize, iiv, insertend, jjv, max_nz, nonzero, release_nodestuff(), srt_hash, and xadj.

Referenced by matvec::Model::vce_emreml(), matvec::Model::vce_gibbs(), and ~SparseMatrix(). 

00504 {
00505    if (adjncy) {delete [] adjncy; adjncy = 0;}
00506    if (xadj) {delete [] xadj; xadj = 0;}
00507    if (iiv) {delete [] iiv; iiv = 0;}
00508    if (jjv) {delete [] jjv; jjv = 0;}
00509    if (aav) {delete [] aav; aav = 0;}
00510    if (hash_srt) {delete [] hash_srt; hash_srt = 0;}
00511    if (srt_hash) {delete [] srt_hash; srt_hash = 0;}
00512    release_nodestuff();
00513    dim=0;
00514    hsize=0;
00515    nonzero=0;
00516    max_nz=0;
00517    insertend=0;
00518 }

void matvec::SparseMatrix::release_nodestuff (  void    ) 
 

Definition at line 520 of file sparsematrix.cpp.

References matvec::Vector< double >::clear(), matvec::Vector< unsigned >::clear(), matvec::Vector< NodeStruct >::clear(), diag, elim_done, factor_done, initial_lij_done, initialize_node_list_done, inv_order, mask, node_list, and order.

Referenced by release(), and resize(). 

00521 {
00522    node_list.clear();
00523    mask.clear();
00524    order.clear();
00525    inv_order.clear();
00526    diag.clear();
00527    initialize_node_list_done=0;
00528    elim_done=0;
00529    factor_done=0;
00530    initial_lij_done=0;
00531 }

void matvec::SparseMatrix::reorder (  const int    path = 2 ) 
 

Re-order elements.

Parameters:
path  an integer path specification; values and their meanings are -
  • 1 find minimum degree ordering only
  • 2 find minimum degree ordering and reorder symmetrically stored matrix (used when only the nonzero entries in the upper triangle of m are being stored)
  • 3 reorder symmetrically stored matrix as specified by input permutation (used when an ordering has already been determined and only the nonzero entries in the upper triangle of m are being stored)
  • 4 same as 2 but put diagonal entries at start of each row
  • 5 same as 3 but put diagonal entries at start of each row

Definition at line 833 of file sparsematrix.cpp.

References initial_lij(), initialize_node_list(), and minfilsymelm_node_list().

Referenced by matvec::Model::graph_log_likelihood_peeling(), matvec::Model::log_likelihood_peeling(), and matvec::Model::vce_emreml_multi_trait(). 

00834 {
00835    if (path ==1) {
00836      initialize_node_list();
00837      minfilsymelm_node_list();
00838    }
00839    else if (path==2) {
00840      initialize_node_list();
00841      minfilsymelm_node_list();
00842      initial_lij();
00843    }
00844    else if (path==3) {
00845      initial_lij();
00846    }
00847    else if (path==4) {
00848      initialize_node_list();
00849      minfilsymelm_node_list();
00850      initial_lij();
00851   }
00852    else if (path==5) {
00853      initial_lij();
00854    }
00855 }

int matvec::SparseMatrix::reset (  const int    ii,
const int    jj,
const double    value
) 
 

Reset the element (i,j) to the value.

Definition at line 1018 of file sparsematrix.cpp.

References a(), aav, dim, matvec::Session::epsilon, hsize, ia(), iiv, ja(), jjv, and matvec::SESSION. 

01019 {
01020   int i=ii;
01021   int j=jj;
01022 
01023    if ( dim == 0 ) throw exception("SparseMatrix::reset(): SparseMatrix is empty");
01024    if (i<=0 || j<=0 ) throw exception("SparseMatrix::reset(i,j,a), range error");
01025    if ( i>j){
01026      i=jj;
01027      j=ii;
01028      //throw exception("SparseMatrix::reset(i,j,a): i<j only upper triangular allowed");
01029    }
01030    
01031    if (fabs(value) <= SESSION.epsilon && i != j) return(1);
01032    int *ia   = iiv-1;
01033    int *ja   = jjv-1;
01034    double *a = aav-1;
01035    unsigned new_ill, ill;
01036    ill =  1 + (433*i+53*j) % hsize;
01037    for (unsigned k=1; k<=200; k++) {
01038       if (ia[ill] == 0) {       // first time store
01039          ia[ill] = i;
01040          ja[ill] = j;
01041          a[ill] = value;
01042          return(1);
01043       }
01044       else if (ia[ill] == i && ja[ill] == j) {
01045          a[ill] = value;      // the previous value has been reset
01046          return(1);
01047       }
01048       else {
01049          new_ill = (ill+640) % hsize + 1;
01050          ill = (new_ill == ill) ? (ill+1)% hsize + 1 : new_ill;
01051       }
01052    }
01053    throw exception("SparseMatrix::insert(): max_nz too small");
01054    return(0);
01055 }

SparseMatrix & matvec::SparseMatrix::resize (  const int    n,
const int    maxnz
) 
 

Resize the object with n by n and maximum non-zero elements maxnz.

if n and maxnz are the same as the previous ones, then we assume that the matrix has the same pattern of non-zeros. So, we will just change the numerical values using initial_lij and then factor etc. This is useful for DFREML

Definition at line 288 of file sparsematrix.cpp.

References aav, dim, factor_done, hash_srt, hsize, iiv, initial_lij_done, insert(), insertend, jjv, logdet_done, max_nz, matvec::next_prime(), nonzero, release_nodestuff(), solver_name, and srt_hash.

Referenced by matvec::GLMM::ainv(), matvec::GLMM::build_hInv(), copyfrom(), matvec::Model::copyfrom(), matvec::Model::fitdata(), matvec::Model::genotype_dist_peeling(), matvec::Model::graph_log_likelihood_peeling(), matvec::Model::log_likelihood_peeling(), matvec::Model::release_mme(), matvec::Population::setup_m_ww(), matvec::Model::setup_mme(), matvec::GLMM::setup_mme(), matvec::Model::setup_ww_single_trait(), matvec::Model::vce_emreml_multi_trait(), matvec::Model::vce_emreml_single_trait(), and matvec::Model::vce_gibbs(). 

00289 {
00290   // std::cout << "\n n " << n << " dim " << dim << " max_nz " << max_nz << " maxnz " << maxnz <<"\n";
00291    if ( dim != n || max_nz != maxnz) {
00292       dim     = n;
00293       max_nz  = maxnz;
00294       if (max_nz == 0) {
00295          hsize = 0;
00296       }
00297       else {
00298          hsize = static_cast<unsigned>(next_prime(static_cast<long>(1.25*max_nz + 50)));
00299       }
00300       if (iiv) {
00301         delete [] iiv;
00302         iiv=0;
00303       }
00304       if (jjv) {
00305         delete [] jjv;
00306         jjv=0;
00307       }
00308       if (aav) {
00309         delete [] aav;
00310         aav=0;
00311       }
00312       if (hash_srt) {
00313         delete [] hash_srt;
00314         hash_srt=0;
00315       }
00316       if (srt_hash) {
00317         delete [] srt_hash;
00318         srt_hash=0;
00319       }
00320 
00321       iiv = new int [hsize];
00322       jjv = new int [hsize];
00323       aav = new double [hsize];
00324       hash_srt = new int [hsize+1];
00325       srt_hash = new int [hsize+1];
00326       nonzero=0;
00327       // std::cout <<"\n release nodes \n";
00328       release_nodestuff(); 
00329       solver_name = 'N';
00330    }  
00331    memset(iiv,'\0',sizeof(int)*hsize);
00332    memset(jjv,'\0',sizeof(int)*hsize); 
00333    if(hsize){
00334      memset(hash_srt,'\0',sizeof(int)*(hsize+1));
00335      memset(srt_hash,'\0',sizeof(int)*(hsize+1));
00336    }
00337    memset(aav,'\0',sizeof(double)*hsize);
00338    initial_lij_done=0;
00339    factor_done=0;
00340    logdet_done = 0;
00341    insertend   = 0;
00342    nonzero     = 0;
00343 
00344    for (unsigned i=1; i<=dim; i++) insert(i,i,0.0);
00345    return *this;
00346 }

Vector< double > matvec::SparseMatrix::row (  const int    ithrow ) 
 

Retrieve the i'th row. Note that the first row is 0.

Definition at line 379 of file sparsematrix.cpp.

References a(), aav, matvec::Vector< T >::begin(), dim, ia(), iiv, ja(), and jjv.

Referenced by initial_lij(), initialize_node_list(), and solvrow(). 

00380 {
00381    if ( dim == 0 ) throw exception("SparseMatrix is empty");
00382    if (ithrow<1 || ithrow>dim) throw exception("SparseMatrix.row(k): out of range");
00383    Vector<double> out(dim);
00384    int *ia = iiv-1;
00385    int *ja = jjv-1;
00386    double *a = aav-1;
00387    double *o_pt = out.begin()-1;
00388    for (unsigned j=ia[ithrow]; j<ia[ithrow+1]; j++) o_pt[ja[j]] = a[j];
00389    return out;
00390 }

void matvec::SparseMatrix::save (  const std::string &    fname,
const int    io_mode = std::ios::out
)  const
 

Save the contents into a disk file.

Definition at line 1151 of file sparsematrix.cpp.

References aav, dim, iiv, and jjv. 

01152 {
01153 
01154    std::ofstream ofs;
01155    ofs.open(fname.c_str(),(OpenModeType)io_mode);
01156    if (!ofs) throw exception(std::string("SparseMatrix::save(): cannot open file:") + fname);
01157 
01158    ofs.precision(16);
01159    ofs.setf(std::ios::unitbuf);
01160    int i,j,ii,jj,jjj,last_j;
01161    for (i=0; i<dim; i++) {
01162       ii = i + 1; last_j = iiv[ii] - 1;
01163       for (j = iiv[i]; j<=last_j; j++) {
01164          jjj = j-1;
01165          jj = jjv[jjj];
01166          ofs << ii << " " << jj << " " << aav[jjj] << "\n";
01167          if (jj != ii) ofs << jj << " " << ii << " " << aav[jjj] << "\n";
01168 
01169       }
01170    }
01171    ofs.close();
01172    ofs.unsetf(std::ios::unitbuf);
01173 }

Vector< double > matvec::SparseMatrix::solv (  Vector< double >    rhs ) 
 

Definition at line 1794 of file sparsematrix.cpp.

References matvec::NodeStruct::a, matvec::NodeStruct::adj_list, matvec::Vector< NodeStruct >::begin(), matvec::Vector< unsigned >::begin(), matvec::Vector< T >::begin(), matvec::Vector< double >::begin(), diag, dim, factor_done, node_list, order, matvec::Vector< T >::reserve(), and matvec::Vector< unsigned >::size(). 

01795 {
01796   if (!factor_done) throw exception("SparseMatrix::factor must be called first");
01797   unsigned i,j;
01798   unsigned coli,rowj,corowj;
01799   unsigned degi;
01800   double diagi,soli;
01801   Vector<double> sol,fsol;
01802   unsigned *prows;
01803   double *pa;
01804   NodeStruct *nodes, *node;
01805 
01806   sol.reserve(dim);
01807   fsol.reserve(dim);
01808 
01809   double   *pdiag  = diag.begin();  pdiag--;
01810   double   *psol   = sol.begin();    psol--;
01811   double   *pfsol  = fsol.begin();   pfsol--;
01812   double   *prhs   = rhs.begin();    prhs--;
01813   unsigned *porder = order.begin(); porder--;
01814 
01815 
01816   // forward elimination
01817   // loop for getting solution 1..dim
01818   nodes = node_list.begin(); nodes--;
01819   for (i=1;i<=dim;i++) {
01820     coli  = porder[i];
01821     diagi = pdiag[coli];
01822     node  = &nodes[coli];
01823     prows = node->adj_list.begin();
01824     pa    = node->a.begin();
01825 
01826     if (diagi == 0) { //lineraly dependent equation
01827       psol[coli] = 0.0;
01828     }
01829     else {
01830       soli = psol[coli] = prhs[coli]/diagi;
01831       // adjust rhs for this column
01832       // note: that diag and rhs are in original order
01833       //       but, the row subscripts of L are  in mindeg order
01834 
01835 
01836       degi = node->adj_list.size();
01837       for (j=0;j<degi;j++) {
01838         rowj = porder[prows[j]];
01839         prhs[rowj] -= pa[j]*soli;
01840       }
01841     }
01842   }
01843   // backward substitution
01844   // loop for getting solution dim..1
01845   for (i=dim;i>=1;i--) {
01846     coli = porder[i];
01847     diagi = pdiag[coli];
01848     node  = &nodes[coli];
01849     prows = node->adj_list.begin();
01850     pa    = node->a.begin();
01851     if (diagi == 0) { //lineraly dependent equation
01852       pfsol[coli] = 0.0;
01853     }
01854     else {
01855       // adjust this rhs for all previous solutions
01856       degi = node->adj_list.size();  
01857       soli=psol[coli];
01858       for (j=0;j<degi;j++) {
01859         rowj = porder[prows[j]];
01860         soli -=  pa[j]*pfsol[rowj];
01861       }
01862       pfsol[coli] = soli/diagi;
01863     }
01864   }
01865   return fsol;
01866 }

int matvec::SparseMatrix::solv (  double *    x,
const double *    rhs
) 
 

Definition at line 1868 of file sparsematrix.cpp.

References matvec::NodeStruct::a, matvec::NodeStruct::adj_list, matvec::Vector< NodeStruct >::begin(), matvec::Vector< unsigned >::begin(), matvec::Vector< double >::begin(), matvec::Vector< T >::begin(), diag, dim, factor_done, node_list, order, matvec::Vector< T >::reserve(), and matvec::Vector< unsigned >::size().

Referenced by solve(). 

01869 {
01870   if (!factor_done) throw exception("SparseMatrix::factor must be called first");
01871   unsigned i,j;
01872   unsigned coli,rowj;
01873   unsigned degi;
01874   double diagi,soli;
01875   Vector<double> sol,temp_rhs;
01876   unsigned *prows;
01877   double *pa;
01878   NodeStruct *nodes, *node;
01879 
01880   sol.reserve(dim);
01881   temp_rhs.reserve(dim);
01882   memcpy(temp_rhs.begin(),rhs,sizeof(double)*dim);
01883 
01884   double   *pdiag  = diag.begin();     pdiag--;
01885   double   *psol   = sol.begin();       psol--;
01886   unsigned *porder = order.begin();    porder--;
01887   double   *prhs   = temp_rhs.begin();  prhs--;
01888   double   *pfsol; 
01889   pfsol = x; pfsol--;
01890  
01891 
01892   
01893   // forward elimination
01894   // loop for getting solution 1..dim
01895   nodes = node_list.begin(); nodes--;
01896   for (i=1;i<=dim;i++) {
01897     coli  = porder[i];
01898     diagi = pdiag[coli];
01899     node  = &nodes[coli];
01900     prows = node->adj_list.begin();
01901     pa    = node->a.begin();
01902 
01903     if (diagi == 0) { //lineraly dependent equation
01904       psol[coli] = 0.0;
01905     }
01906     else {
01907       soli = psol[coli] = prhs[coli]/diagi;
01908       // adjust rhs for this column
01909       // note: that diag and rhs are in original order
01910       //       but, the row subscripts of L are  in mindeg order
01911 
01912 
01913       degi = node->adj_list.size();
01914       for (j=0;j<degi;j++) {
01915         rowj = porder[prows[j]];
01916         prhs[rowj] = prhs[rowj] - pa[j]*soli;
01917       }
01918     }
01919   }
01920   // backward substitution
01921   // loop for getting solution dim..1
01922   for (i=dim;i>=1;i--) {
01923     coli = porder[i];
01924     diagi = pdiag[coli];
01925     node  = &nodes[coli];
01926     prows = node->adj_list.begin();
01927     pa    = node->a.begin();
01928     if (diagi == 0) { //lineraly dependent equation
01929       pfsol[coli] = 0.0;
01930     }
01931     else {
01932       // adjust this rhs for all previous solutions
01933       degi = node->adj_list.size();
01934       for (j=0;j<degi;j++) {
01935         rowj = porder[prows[j]];
01936         psol[coli] = psol[coli] - pa[j]*pfsol[rowj];
01937       }
01938       pfsol[coli] = psol[coli]/diagi;
01939     }
01940   }
01941   return 1;
01942 }

int matvec::SparseMatrix::solve (  Vector< double > &    x,
const Vector< double > &    b,
const std::string &    method = "ysmp",
double    relax = 1.0,
double    stopval = 0.001,
int    mxiter = 1000
) 
 

Obtain a solution with the right-hand side b.

Parameters:
x  a solution
b  right hand side
method  solver
relax  relaxation factor
stopval  criteria for stop iteration
mxiter  maximum iteration allowed

Definition at line 752 of file sparsematrix.cpp.

References aav, matvec::Vector< T >::begin(), close(), dim, matvec::Session::epsilon, factor(), factor_done, matvec::gs_ioc(), iiv, jjv, matvec::Vector< T >::reserve(), matvec::SESSION, matvec::Vector< T >::size(), and solv(). 

00754 {
00755    if (dim != b.size()) throw exception("SparseMatrix::solve(): size not conformable");
00756    if (x.size() < dim ) x.reserve(dim);
00757    if (method == "ysmp" || method == "ysmp1") {
00758       if (!factor_done && !factor()) throw exception("SparseMatrix::solve- returned with error from factor");
00759       x = solv(b);
00760       return 1;
00761    }
00762    else if (method == "ioc") {
00763       close();
00764       unsigned ir=0;
00765       memset(x.begin(),'\0',sizeof(double)*dim);
00766       gs_ioc(iiv-1,jjv-1,aav-1,b.begin()-1,x.begin()-1,dim,relax,stopval,mxiter,SESSION.epsilon,ir);
00767       if (!ir) throw exception("SparseMatrix::solve(): not converged");
00768       return 1;
00769    }
00770    else if (method == "iod") {
00771       throw exception("SparseMatrix::solve(): not available yet");
00772    }
00773    else {
00774       throw exception("SparseMatrix::solve(): no such solver");
00775   }
00776   return 0;
00777 }

int matvec::SparseMatrix::solve (  double *    x,
const double *    b,
const std::string &    method = "ysmp",
double    relax = 1.0,
double    stopval = 0.001,
int    mxiter = 1000
) 
 

Obtain a solution with the right-hand side b.

Parameters:
x  a solution
b  right hand side
method  solver
relax  relaxation factor
stopval  criteria for stop iteration
mxiter  maximum iteration allowed
Warning:
x and b must have at least this->dim elements x[0] and b[0] are the first elements.

Definition at line 792 of file sparsematrix.cpp.

References aav, close(), dim, matvec::Session::epsilon, factor(), factor_done, matvec::gs_ioc(), iiv, jjv, matvec::SESSION, solv(), and matvec::warning().

Referenced by matvec::Model::blup(), matvec::GLMM::contrast(), matvec::Model::CovMat(), matvec::Model::genotype_dist_peeling(), matvec::Model::genotypic_value_peeling(), matvec::Model::get_blupsol(), matvec::GLMM::glim(), matvec::Model::graph_log_likelihood_peeling(), matvec::GLMM::log_likelihood(), matvec::Model::log_likelihood_peeling(), matvec::GLMM::residual(), matvec::Model::restricted_log_likelihood(), matvec::GLMM::restricted_log_likelihood(), matvec::GLMM::SSQCP(), matvec::Model::uAu_trCA(), matvec::Model::vce_emreml_multi_trait(), and matvec::Model::vce_emreml_single_trait(). 

00794 {
00795    if (method == "ysmp" || method == "ysmp1") {
00796       if (!factor_done && !factor()) throw exception("SparseMatrix::solve() returned with error from factor");
00797       solv(x,b);
00798       return 1;
00799     }
00800     else if (method == "ioc") {
00801       close();
00802       unsigned ir=0;
00803       memset(x,'\0',sizeof(double)*dim);
00804       gs_ioc(iiv-1,jjv-1,aav-1,b-1,x-1,dim,relax,stopval,
00805                mxiter,SESSION.epsilon,ir);
00806       if (!ir) warning("SparseMatrix::solve(): not converged");
00807       return 1;
00808     }
00809     else if (method == "iod") {
00810       throw exception("SparseMatrix::solve(): solver not available yet");
00811    }
00812    else {
00813       throw exception("SparseMatrix::solve():no such solver");
00814    }
00815 }

int matvec::SparseMatrix::solvrow (  double *    x,
const double *    rhs,
unsigned    row,
Vector< double > &    sol,
Vector< double > &    temp_rhs,
unsigned *    row_pt,
unsigned *    col_pt
) 
 

Definition at line 2012 of file sparsematrix.cpp.

References matvec::NodeStruct::a, matvec::NodeStruct::adj_list, matvec::Vector< NodeStruct >::begin(), matvec::Vector< unsigned >::begin(), matvec::Vector< T >::begin(), matvec::Vector< double >::begin(), diag, dim, factor_done, inv_order, node_list, order, row(), and matvec::Vector< unsigned >::size().

Referenced by matvec::GLMM::build_hInv(). 

02012                                                                                                                                                   {
02013   unsigned i,j,allzero;
02014   unsigned coli,rowj,corowj;
02015   unsigned degi;
02016   double diagi,soli;
02017   //Vector<double> sol,temp_rhs;
02018   unsigned *prows;
02019   double *pa;
02020   NodeStruct *nodes, *node;
02021   
02022   
02023   if (!factor_done) {
02024     throw exception("SparseMatrix::factor must be called first");
02025     return 0;
02026   }
02027   //  sol.reserve(dim);
02028   //  temp_rhs.reserve(dim);
02029   // unsigned *row_pt=new unsigned [dim];
02030   unsigned *row_mask=row_pt;
02031   //  unsigned *col_pt=new unsigned [dim];
02032   unsigned *col_mask=col_pt;
02033   row_mask--;
02034   memset(row_pt,'\0',dim*sizeof(unsigned));
02035   col_mask--;
02036   memset(col_pt,'\0',dim*sizeof(unsigned));
02037   allzero=1;
02038   double   *pdiag  = diag.begin();     pdiag--;
02039   double   *psol   = sol.begin();       psol--;
02040   unsigned *porder = order.begin();    porder--;
02041   unsigned *pinv_order = inv_order.begin();    pinv_order--;
02042   double   *prhs   = temp_rhs.begin();  prhs--;
02043   double   *pfsol; 
02044   pfsol = x; pfsol--;
02045   memset(temp_rhs.begin(),'\0',dim*sizeof(double));
02046   //  prhs[porder[row]]=1;
02047   prhs[row]=1;
02048   //  memcpy(temp_rhs.begin(),rhs,dim*sizeof(double));
02049   memset(x,'\0',dim*sizeof(double));
02050 
02051   
02052   // forward elimination
02053   // loop for getting solution 1..dim
02054   nodes = node_list.begin(); nodes--;
02055   row_mask[row]=1;
02056   col_mask[row]=1;
02057   for (i=row;i<=dim;i++) {
02058     if(row_mask[i]) {
02059       coli  = porder[i];
02060       diagi = pdiag[coli];
02061       node  = &nodes[coli];
02062       prows = node->adj_list.begin();
02063       pa    = node->a.begin(); 
02064       //row_mask[i]=1;
02065       if (diagi != 0)  {
02066         //      soli = psol[coli] = prhs[coli]/diagi;
02067         soli = psol[i] = prhs[i]/diagi;
02068         row_mask[i]=1;
02069         // adjust rhs for this column
02070         // note: that diag and rhs are in original order
02071         //       but, the row subscripts of L are  in mindeg order
02072         
02073         
02074         degi = node->adj_list.size();
02075         for (j=0;j<degi;j++) {
02076           row_mask[prows[j]]=1;
02077           //      if(prows[j] == row) col_mask[i]=1;
02078           //      rowj = porder[prows[j]];
02079           //      prhs[rowj] -=  pa[j]*soli;
02080           prhs[prows[j]] -=  pa[j]*soli;
02081         }
02082       }
02083     }
02084   }
02085   // backward substitution
02086   // loop for getting solution dim..1
02087   for (i=dim;i>=row;i--) {
02088     if(row_mask[i]) {
02089       coli = porder[i];
02090       diagi = pdiag[coli];
02091       if (diagi != 0 )  {
02092         node  = &nodes[coli];
02093         prows = node->adj_list.begin();
02094         pa    = node->a.begin();
02095         // adjust this rhs for all previous solutions
02096         degi = node->adj_list.size();
02097         //      soli=psol[coli];
02098 #pragma ivdep
02099         for (j=0;j<degi;j++) {
02100           //      rowj = porder[prows[j]];
02101           //      psol[coli] -= pa[j]*pfsol[rowj];
02102           psol[i] -= pa[j]*pfsol[prows[j]];
02103         }
02104         //      pfsol[coli] = psol[coli]/diagi;
02105         pfsol[i] = psol[i]/diagi;
02106       }
02107     }
02108   }
02109   //  for(i=1;i<=dim;i++){
02110   //    psol[porder[i]]=pfsol[i];
02111   //  }
02112   //  delete [] row_pt;
02113   //  delete [] col_pt;
02114   return 1;
02115 }

Friends And Related Function Documentation

friend class GLMM [friend]
 

Definition at line 55 of file sparsematrix.h.

Vector<double> operator * (  const SparseMatrix &    A,
const Vector< double > &    V
)  [friend]
 

SparseMatrix times vector operation.

Definition at line 395 of file sparsematrix.cpp. 

00396 {
00397   if ( A.dim == 0) throw exception("SparseMatrix is empty");
00398   if (A.dim != v.size()) throw exception("SparseMatrix*Vector<double>, size not conformable");
00399   Vector<double> vec_out;
00400   int n = A.dim;
00401   int *A_ia = A.iiv-1;   /* offset by 1 */
00402   int *A_ja = A.jjv-1;
00403   double *A_a = A.aav-1;
00404   double *v_ve = v.begin()-1;
00405   double tmp,x;
00406   unsigned k,i,j;
00407 
00408   vec_out.reserve(n);
00409   double *out = vec_out.begin();
00410 
00411   out--; // offset by 1
00412   int ksrt=0;       //SDK
00413   for(ksrt=0;ksrt<A.nonzero;ksrt++) { //SDK
00414     k=A.srt_hash[ksrt];  //SDK
00415   //  for (k=1; k<=A.hsize; k++) {  //SDK
00416     if (i = A_ia[k]) {
00417         j = A_ja[k];
00418         x = A_a[k];
00419         out[i] = out[i] + x*v_ve[j];
00420         if (j>i) {
00421           out[j] = out[j] + x*v_ve[i];
00422         }
00423     }
00424   }
00425   return vec_out;
00426 }

SparseMatrix operator * (  const double    s,
const SparseMatrix &    A
)  [friend]
 

Definition at line 499 of file sparsematrix.cpp. 

00500 {return operator*(A,s);}

SparseMatrix operator * (  const SparseMatrix &    A,
const double    s
)  [friend]
 

Multiply SparseMatrix object by s.

Definition at line 475 of file sparsematrix.cpp. 

00476 {
00477    if ( A.dim == 0 ) throw exception("SparseMatrix is empty");
00478 
00479    unsigned new_dim = A.dim;
00480    unsigned new_nz = A.nonzero;
00481    SparseMatrix B(new_dim,A.max_nz);
00482    B.nonzero = new_nz;
00483    memcpy(B.iiv,A.iiv,sizeof(int)*A.hsize);
00484    memcpy(B.jjv,A.jjv,sizeof(int)*A.hsize);
00485    if(A.hsize){
00486      memcpy(B.hash_srt,A.hash_srt,sizeof(int)*(A.hsize+1));
00487      memcpy(B.srt_hash,A.srt_hash,sizeof(int)*(A.hsize+1));
00488    }
00489    double *A_pt = A.aav;
00490    double *New_pt = B.aav;
00491    for (unsigned i=0; i<A.hsize; i++) {
00492      if (A.iiv[i] > 0) {
00493        New_pt[i] = A_pt[i]*s;
00494      }
00495    }
00496    return B;
00497 }

int operator!= (  const SparseMatrix &    A,
const SparseMatrix &    B
)  [friend]
 

Not equal operator.

Definition at line 373 of file sparsematrix.cpp. 

00374 {return !operator==(A,B);}

std::ostream& operator<< (  std::ostream &    stream,
const SparseMatrix &    A
)  [friend]
 

int operator== (  const SparseMatrix &    A,
const SparseMatrix &    B
)  [friend]
 

Relational operator. Return 1 if A and B are identical, otherwise return 0.

Definition at line 352 of file sparsematrix.cpp. 

00353 {
00354    if (A.dim != B.dim || A.nonzero != B.nonzero) return 0;
00355    unsigned n = A.dim;
00356    unsigned m = A.nonzero;
00357    if (m*n == 0) return 1;
00358 
00359    unsigned i;
00360    double *A_a = A.aav-1;
00361    double *B_a = B.aav-1;
00362 
00363    // one loop instead of the two loops of Tianlin
00364    for (i=1; i<=m; i++) {
00365      if (A_a[i] != B_a[i])  return 0;
00366    }
00367    return 1;
00368 }

Member Data Documentation

double* matvec::SparseMatrix::aav [protected]
 

Definition at line 65 of file sparsematrix.h.

Referenced by a(), add_diag(), matvec::GLMM::build_hInv(), close(), copyfrom(), getaij(), gibbs_iterate(), initialize_node_list(), insert(), mv(), matvec::operator *(), matvec::operator==(), q(), qTLW(), release(), reset(), resize(), row(), save(), solve(), and SparseMatrix().

idxtype * matvec::SparseMatrix::adjncy [protected]
 

Definition at line 64 of file sparsematrix.h.

Referenced by minfilsymelm_node_list(), release(), and SparseMatrix().

Vector<double> matvec::SparseMatrix::diag [protected]
 

Definition at line 61 of file sparsematrix.h.

Referenced by factor(), initial_lij(), logdet(), minfilsymelm_node_list(), release_nodestuff(), solv(), and solvrow().

unsigned matvec::SparseMatrix::dim [protected]
 

Definition at line 62 of file sparsematrix.h.

Referenced by matvec::GLMM::build_hInv(), close(), copyfrom(), dense(), display(), display_inv_order(), display_lij(), display_node_list(), display_order(), factor(), getaij(), gibbs_iterate(), initial_lij(), initialize_node_list(), insert(), logdet(), minfilsymelm_node_list(), mv(), num_rows(), matvec::operator *(), operator()(), matvec::operator==(), q(), qTLW(), release(), reset(), resize(), row(), save(), solv(), solve(), solvrow(), and SparseMatrix().

unsigned matvec::SparseMatrix::elim_done [protected]
 

Definition at line 69 of file sparsematrix.h.

Referenced by copyfrom(), initial_lij(), minfilsymelm_node_list(), release_nodestuff(), and SparseMatrix().

unsigned matvec::SparseMatrix::factor_done [protected]
 

Definition at line 69 of file sparsematrix.h.

Referenced by matvec::GLMM::build_hInv(), copyfrom(), factor(), logdet(), release_nodestuff(), resize(), solv(), solve(), solvrow(), and SparseMatrix().

int * matvec::SparseMatrix::hash_srt [protected]
 

Definition at line 66 of file sparsematrix.h.

Referenced by copyfrom(), insert(), matvec::operator *(), release(), resize(), and SparseMatrix().

unsigned matvec::SparseMatrix::hsize [protected]
 

Definition at line 62 of file sparsematrix.h.

Referenced by matvec::GLMM::build_hInv(), close(), copyfrom(), getaij(), insert(), matvec::operator *(), release(), reset(), resize(), and SparseMatrix().

int* matvec::SparseMatrix::iiv [protected]
 

Definition at line 63 of file sparsematrix.h.

Referenced by add_diag(), matvec::GLMM::build_hInv(), close(), copyfrom(), getaij(), gibbs_iterate(), ia(), initialize_node_list(), insert(), mv(), matvec::operator *(), q(), qTLW(), release(), reset(), resize(), row(), save(), solve(), and SparseMatrix().

Vector<double> matvec::SparseMatrix::info_vec
 

Definition at line 75 of file sparsematrix.h.

Referenced by copyfrom(), factor(), matvec::Model::graph_log_likelihood_peeling(), matvec::Model::log_likelihood_peeling(), matvec::Model::restricted_log_likelihood(), matvec::GLMM::restricted_log_likelihood(), and SparseMatrix().

unsigned matvec::SparseMatrix::initial_lij_done [protected]
 

Definition at line 69 of file sparsematrix.h.

Referenced by copyfrom(), factor(), initial_lij(), release_nodestuff(), resize(), and SparseMatrix().

unsigned matvec::SparseMatrix::initialize_node_list_done [protected]
 

Definition at line 69 of file sparsematrix.h.

Referenced by copyfrom(), initialize_node_list(), minfilsymelm_node_list(), release_nodestuff(), and SparseMatrix().

unsigned matvec::SparseMatrix::insertend [protected]
 

Definition at line 62 of file sparsematrix.h.

Referenced by close(), copyfrom(), qTLW(), release(), resize(), and SparseMatrix().

Vector<unsigned> matvec::SparseMatrix::inv_order [protected]
 

Definition at line 60 of file sparsematrix.h.

Referenced by matvec::GLMM::build_hInv(), display_inv_order(), minfilsymelm_node_list(), release_nodestuff(), and solvrow().

int * matvec::SparseMatrix::iperm [protected]
 

Definition at line 63 of file sparsematrix.h.

int * matvec::SparseMatrix::jjv [protected]
 

Definition at line 63 of file sparsematrix.h.

Referenced by add_diag(), matvec::GLMM::build_hInv(), close(), copyfrom(), getaij(), gibbs_iterate(), initialize_node_list(), insert(), ja(), mv(), matvec::operator *(), q(), qTLW(), release(), reset(), resize(), row(), save(), solve(), and SparseMatrix().

unsigned matvec::SparseMatrix::logdet_done [protected]
 

Definition at line 70 of file sparsematrix.h.

Referenced by irank(), logdet(), resize(), and SparseMatrix().

double matvec::SparseMatrix::logdeterm
 

Definition at line 74 of file sparsematrix.h.

Referenced by logdet().

Vector<unsigned> matvec::SparseMatrix::mask [protected]
 

Definition at line 58 of file sparsematrix.h.

Referenced by elim_node(), initialize_node_list(), minfilsymelm_node_list(), and release_nodestuff().

Vector<unsigned> matvec::SparseMatrix::masklist [protected]
 

Definition at line 58 of file sparsematrix.h.

Referenced by elim_node(), and minfilsymelm_node_list().

unsigned matvec::SparseMatrix::max_nz [protected]
 

Definition at line 62 of file sparsematrix.h.

Referenced by copyfrom(), matvec::operator *(), release(), resize(), and SparseMatrix().

Vector<NodeStruct> matvec::SparseMatrix::node_list [protected]
 

Definition at line 57 of file sparsematrix.h.

Referenced by display_lij(), display_node_list(), elim_node(), factor(), initial_lij(), initialize_node_list(), minfilsymelm_node_list(), release_nodestuff(), solv(), and solvrow().

unsigned matvec::SparseMatrix::nonzero [protected]
 

Definition at line 62 of file sparsematrix.h.

Referenced by close(), copyfrom(), initialize_node_list(), insert(), minfilsymelm_node_list(), mv(), nz(), matvec::operator *(), matvec::operator==(), q(), release(), resize(), and SparseMatrix().

unsigned matvec::SparseMatrix::nsp [protected]
 

Definition at line 62 of file sparsematrix.h.

Vector<unsigned> matvec::SparseMatrix::order [protected]
 

Definition at line 59 of file sparsematrix.h.

Referenced by matvec::GLMM::build_hInv(), display_node_list(), display_order(), factor(), initial_lij(), minfilsymelm_node_list(), release_nodestuff(), solv(), and solvrow().

int * matvec::SparseMatrix::perm [protected]
 

Definition at line 63 of file sparsematrix.h.

unsigned matvec::SparseMatrix::rank [protected]
 

Definition at line 71 of file sparsematrix.h.

Referenced by irank(), and logdet().

double * matvec::SparseMatrix::rsp [protected]
 

Definition at line 65 of file sparsematrix.h.

char matvec::SparseMatrix::solver_name [protected]
 

Definition at line 67 of file sparsematrix.h.

Referenced by copyfrom(), resize(), and SparseMatrix().

int* matvec::SparseMatrix::srt_hash [protected]
 

Definition at line 66 of file sparsematrix.h.

Referenced by copyfrom(), initialize_node_list(), insert(), mv(), matvec::operator *(), q(), release(), resize(), and SparseMatrix().

idxtype* matvec::SparseMatrix::xadj [protected]
 

Definition at line 64 of file sparsematrix.h.

Referenced by minfilsymelm_node_list(), release(), and SparseMatrix().

The documentation for this class was generated from the following files:
Generated on Thu Jun 16 17:14:50 2005 for Matvec by doxygen1.2.16