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

matvec::SparseBGMatrix Class Reference

#include <sparsebgmatrix.h>

List of all members.

Detailed Description

a SparseBGMatrix is a sparse matrix of BG objects

See also:
Matrix

Definition at line 57 of file sparsebgmatrix.h.

Public Methods

 SparseBGMatrix (void)
 SparseBGMatrix (int n, int max_nz)
 SparseBGMatrix (const SparseBGMatrix &A)
 ~SparseBGMatrix (void)
const SparseBGMatrix & operator= (const SparseBGMatrix &A)
BG operator() (const int i, const int j)
SparseBGMatrix & resize (const int n, const int maxnz)
BG q (const Vector< BG > &v1, const Vector< BG > &v2)
BG q (const BG *x, const BG *y)
BG q (const BG *x, const BG *y, const int i1, const int i2)
BG qTLW (const Vector< BG > &v1, const Vector< BG > &v2)
BG qTLW (const BG *x, const BG *y)
BG qTLW (const BG *x, const BG *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)
BG * a (void)
doubleMatrix dense (unsigned r1, unsigned r2, unsigned c1, unsigned c2) const
doubleMatrix dense (void)
Vector< BG > row (const int ithrow)
BG logdet (unsigned &irank)
BG logdet (unsigned *irank)
BG logdet (void)
unsigned irank (void)
int solve (BG *x, const BG *b, const std::string &method="ysmp", double relax=1.0, double stopval=0.001, int mxiter=1000)
int solve (Vector< BG > &x, const Vector< BG > &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 BG value)
int insert (const int i, const int j, const BG value)
void reorder (const int path=2)
void add_diag (const int ibeg, const int iend, const BG r)
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< BG > &v, Vector< BG > &result)
void mv (const BG *v, const unsigned n, BG *result)
BG 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 (BG *x, const BG *rhs)
int solvrow (BG *x, const BG *rhs, unsigned row, Vector< BG > &sol, Vector< BG > &temp_rhs, unsigned *row_pt, unsigned *col_pt)
Vector< BG > solv (Vector< BG > rhs)

Public Attributes

BG logdeterm
Vector< BG > info_vec

Protected Methods

void copyfrom (const SparseBGMatrix &A)

Protected Attributes

Vector< BGNodeStruct > node_list
Vector< unsigned > mask
Vector< unsigned > masklist
Vector< unsigned > order
Vector< unsigned > inv_order
Vector< BG > 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
BG * aav
BG * 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 SparseBGMatrix &A, const SparseBGMatrix &B)
int operator!= (const SparseBGMatrix &A, const SparseBGMatrix &B)
SparseBGMatrix operator * (const SparseBGMatrix &A, const BG s)
SparseBGMatrix operator * (const BG s, const SparseBGMatrix &A)
Vector< BG > operator * (const SparseBGMatrix &A, const Vector< BG > &V)
std::ostream & operator<< (std::ostream &stream, const SparseBGMatrix &A)

Constructor & Destructor Documentation

matvec::SparseBGMatrix::SparseBGMatrix (  void    ) 
 

A constructor to create an SparseBGMatrix object.

Definition at line 140 of file sparsebgmatrix.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< BG >::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::SparseBGMatrix::SparseBGMatrix (  int    n,
int    maxnz
) 
 

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

Definition at line 195 of file sparsebgmatrix.cpp.

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

00196 {
00197    if (n<0 || maxnz<0) throw exception("SparseBGMatrix(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 BG [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(BG)*hsize);
00222    initial_BGarray(aav,hsize,0.0);
00223    if(hsize){
00224      memset(hash_srt,'\0',sizeof(int)*(hsize+1));
00225      memset(srt_hash,'\0',sizeof(int)*(hsize+1));
00226    }
00227    nonzero=0;
00228    BG temp(0,0.0);
00229    for (unsigned i=1; i<=dim; i++) insert(i,i,temp);
00230 }

matvec::SparseBGMatrix::SparseBGMatrix (  const SparseBGMatrix &    A ) 
 

A constructor to create an identical SparseBGMatrix object to A.

Definition at line 167 of file sparsebgmatrix.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::SparseBGMatrix::~SparseBGMatrix (  void    )  [inline]
 

Definition at line 85 of file sparsebgmatrix.h.

References release(). 

00085 {release();}                      // Destructor

Member Function Documentation

BG* matvec::SparseBGMatrix::a (  void    )  [inline]
 

Definition at line 126 of file sparsebgmatrix.h.

Referenced by add_diag(), close(), getaij(), initialize_node_list(), insert(), minfilsymelm_node_list(), q(), qTLW(), reset(), and row(). 

00126 {return aav-1;};

void matvec::SparseBGMatrix::add_diag (  const int    ibeg,
const int    iend,
const BG    r
) 
 

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

Definition at line 735 of file sparsebgmatrix.cpp.

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

00736 {
00737    int   *ia = iiv-1;
00738    int   *ja = jjv-1;
00739    BG *a = aav-1;
00740    int i,j;
00741    for (i=ibeg; i<=iend; i++) {
00742       for (j=ia[i]; j<ia[i+1]; j++) if (ja[j] == i) a[j] += r;
00743    }
00744 }

unsigned matvec::SparseBGMatrix::close (  void    ) 
 

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

Definition at line 1066 of file sparsebgmatrix.cpp.

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

Referenced by solve(). 

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

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

Get a copy from A

Definition at line 235 of file sparsebgmatrix.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 SparseBGMatrix(). 

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

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

Definition at line 128 of file sparsebgmatrix.h.

References dim.

Referenced by display(). 

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

doubleMatrix matvec::SparseBGMatrix::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 1111 of file sparsebgmatrix.cpp.

References matvec::BG::f. 

01112 {
01113   doubleMatrix out(r2-r1+1,c2-c1+1);
01114   int   *ia = iiv-1;
01115   int   *ja = jjv-1;
01116   BG *a = aav-1;
01117 
01118   unsigned i,j,k;
01119 
01120   for (k=1; k<= hsize; k++) {
01121     if (i = ia[k]) {
01122         j = ja[k];
01123         if (i!=j) {
01124           if (i>=r1 && i<=r2 && j>=c1 && j<=c2 ) {
01125             out(i-r1+1,j-c1+1) = a[k].f;
01126           }
01127           if (j>=r1 && j<=r2 && i>=c1 && i<=c2 ) {
01128             out(j-r1+1,i-c1+1) = a[k].f;
01129           }
01130         }
01131         else {
01132           if (i>=r1 && i<=r2) {
01133            out(i-r1+1,i-r1+1) = a[k].f;
01134           }
01135         }
01136     }
01137   }
01138   return out;
01139 }

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

Definition at line 146 of file sparsebgmatrix.h.

References dim, and display(). 

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

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

Display the content.

Definition at line 1144 of file sparsebgmatrix.cpp.

References dense(), and dim.

Referenced by display(). 

01145 {
01146   if (msg != "") std::cout <<"\n     " << msg << "\n";
01147   if (dim == 0) {
01148     std::cout << "      null sparse matrix     \n\n";
01149     return;
01150   }
01151   std::cout << dense(r1,r2,c1,c2);
01152 }

void matvec::SparseBGMatrix::display_inv_order (  void    ) 
 

Definition at line 1981 of file sparsebgmatrix.cpp.

References dim, and inv_order. 

01981                                            {
01982   // displays the order of diagonal nodes
01983   unsigned i;
01984   std::cout << "\n" << "inv_order" << "\n";
01985   for (i=1;i<=dim;i++) {
01986     std::cout << inv_order(i) << std::endl;
01987   }
01988   std::cout << std::endl;
01989 }

void matvec::SparseBGMatrix::display_lij (  void    ) 
 

Definition at line 1957 of file sparsebgmatrix.cpp.

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

01957                                      {
01958   // displays the lij's
01959   unsigned i,j;
01960   unsigned count;
01961 
01962   for (i=1;i<=dim;i++) {
01963     std::cout <<"\n" <<"node: "<<i<<" ";
01964     count = node_list(i).adj_list.size();
01965     for (j=1;j<=count;j++){
01966       std::cout << node_list(i).a(j).f <<" ";
01967     }
01968   }
01969   std::cout << std::endl;
01970 }

void matvec::SparseBGMatrix::display_node_list (  void    ) 
 

Definition at line 1926 of file sparsebgmatrix.cpp.

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

01926                                            {
01927   // displays the nodes in the adjacent list
01928   // using peeling order numbers
01929   unsigned i,j;
01930   unsigned count;
01931 
01932   for (i=1;i<=dim;i++) {
01933     int orderi = order(i);
01934     std::cout <<"\n" <<"node: " << i <<": ";
01935     
01936     count = node_list(orderi).adj_list.size();
01937     for (j=1;j<=count;j++){
01938       std::cout << node_list(orderi).adj_list(j)<<" ";
01939     }
01940   }
01941   std::cout << std::endl;
01942 
01943 //   unsigned i,j;
01944 //   unsigned count;
01945 //   for (i=1;i<=dim;i++) {
01946 //     std::cout <<"\n" <<"node: "<<i<<" ";
01947 //     count = node_list(i).adj_list.size();
01948 //     for (j=1;j<=count;j++){
01949 //       std::cout << node_list(i).adj_list(j) <<" ";
01950 //     }
01951 //   }
01952 //   std::cout << std::endl;
01953 
01954 
01955 
01956 }

void matvec::SparseBGMatrix::display_order (  void    ) 
 

Definition at line 1971 of file sparsebgmatrix.cpp.

References dim, and order. 

01971                                        {
01972   // displays the order of diagonal nodes
01973   unsigned i;
01974   std::cout << "\n" << "order" << "\n";
01975   for (i=1;i<=dim;i++) {
01976     std::cout << order(i) << std::endl;
01977   }
01978   std::cout << std::endl;
01979 }

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

Definition at line 2093 of file sparsebgmatrix.cpp.

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

Referenced by minfilsymelm_node_list(). 

02094 {
02095   
02096   // updates elimination graph for eliminating elimnode
02097   
02098   unsigned j,k;
02099   unsigned naborj,nabork,nabor_nabor;
02100   unsigned degnode,degnabor;
02101   unsigned count,*adjpt,*nadjpt,*aadjpt;
02102   unsigned *maskptr,nmask;
02103   //  Vector<unsigned> masklist(dim);
02104   BGNodeStruct *elimnode, *nabornode;
02105   elimnode = &node_list(Aorder(i));
02106   maskptr=mask.begin();
02107   // loop to update adjacent lists for neigbors of elimnode
02108   degnode = elimnode->adj_list.size();
02109   //  memset(mask.begin(),'\0',dim*sizeof(unsigned)); // SDK
02110 
02111 
02112  
02113   adjpt=elimnode->adj_list.begin();
02114   for(k=1;k<=degnode;k++) {
02115     naborj = *adjpt++;
02116     mask[naborj-1]=1;
02117     masklist[k-1]=naborj-1;
02118   }
02119  
02120   adjpt=elimnode->adj_list.begin();
02121   for (j=1;j<=degnode;j++) { 
02122  
02123     naborj = *adjpt++;
02124     
02125 
02126     nabornode = &node_list(naborj);
02127     nmask=degnode;
02128     degnabor = nabornode->adj_list.size();
02129     nadjpt = nabornode->adj_list.begin();
02130     for (k=0;k<degnabor;k++){
02131       nabor_nabor = *nadjpt++;
02132       if (nabor_nabor != Aorder(i)) { //elim elimnode
02133         if(!mask[nabor_nabor-1]) {  //SDK
02134           masklist[nmask++]=nabor_nabor-1; //SDK
02135         }                      //SDK
02136       }
02137     }
02138     count=nmask-1; //SDK
02139     if (nabornode->adj_list.size() != count) {
02140       nabornode->adj_list.reserve(count);
02141       if (nabornode->adj_list.size() < count) return 0;
02142     } 
02143     int adj_count=0;
02144     nadjpt=nabornode->adj_list.begin();
02145     for(count=0;count<nmask;count++){                  //SDK
02146       if(masklist[count] != (naborj-1)) {
02147         adj_count++;
02148         *nadjpt++ = masklist[count]+1;  //SDK
02149       }
02150     }
02151     if(count != (adj_count+1)) {
02152       std::cout << count << " " << adj_count << "\n" << std::flush;
02153       throw exception("Elim node: not possible");
02154     }//SDK
02155   }
02156   
02157   adjpt=elimnode->adj_list.begin();
02158   for(k=1;k<=degnode;k++) {
02159     naborj = *adjpt++;
02160     mask[naborj-1]=0;
02161   }
02162   return 1;
02163 }

int matvec::SparseBGMatrix::factor (  void    ) 
 

Definition at line 1641 of file sparsebgmatrix.cpp.

References matvec::BGNodeStruct::a, matvec::BGNodeStruct::adj_list, matvec::Vector< BGNodeStruct >::begin(), matvec::Vector< unsigned >::begin(), matvec::Vector< BG >::begin(), matvec::Vector< T >::begin(), diag, dim, factor_done, info_vec, initial_lij(), initial_lij_done, matvec::log(), node_list, order, matvec::Vector< T >::reserve(), matvec::Vector< T >::resize(), matvec::Vector< unsigned >::size(), and matvec::sqrt().

Referenced by factorization(), logdet(), and solve(). 

01641                               {
01642   unsigned i,j,nextj,k,count,m;
01643   unsigned coli,colj,rowm,rowk;
01644   unsigned degi,degj;
01645   unsigned startj;
01646   Vector<BG> temp;
01647   Vector<unsigned> start, contr_cols;
01648   BG lij,diagi;
01649   unsigned *prows;
01650   BG   *plij;
01651   BGNodeStruct *nodes, *node;
01652 
01653   if (!initial_lij_done) {
01654     if(!initial_lij()) return 0;
01655   }
01656   temp.reserve(dim);
01657   start.reserve(dim);
01658   contr_cols.reserve(dim);
01659 
01660 
01661   BG   *ptemp        = temp.begin();       ptemp--;
01662   BG   *pdiag        = diag.begin();       pdiag--;
01663   unsigned *pstart       = start.begin();      pstart--;
01664   unsigned *porder       = order.begin();      porder--;
01665   unsigned *pcontri_cols = contr_cols.begin(); pcontri_cols--;
01666   BG    mult;
01667   unsigned lastpos;
01668   lastpos=0;
01669 
01670   info_vec[0]=0;
01671   info_vec[1]=0;
01672   //for (i=1;i<=dim;i++) {
01673   //  ptemp[i]        = 0.0;
01674   //  pstart[i]       = 0;
01675   //  pcontri_cols[i] = 0;
01676   // }
01677   memset(contr_cols.begin(),'\0',dim*sizeof(unsigned));
01678   memset(start.begin(),'\0',dim*sizeof(unsigned));
01679   //memset(temp.begin(),'\0',dim*sizeof(BG));
01680   temp.resize(dim,0.0);
01681   // loop for calculating column i of lij
01682   nodes = node_list.begin(); nodes--;
01683   for (i=1;i<=dim;i++){
01684     coli = porder[i];
01685     j = pcontri_cols[i];
01686     count = 0;
01687     // loop for getting contributions from columns contributing to column i
01688     // contri_cols is used to point to these columns
01689     while (j) {
01690       startj = pstart[j];
01691       colj  = porder[j];
01692       node  = &nodes[colj];
01693       degj  = node->adj_list.size();
01694       nextj = pcontri_cols[j];
01695       prows = node->adj_list.begin();
01696       plij  = node->a.begin();
01697       lij = plij[startj];
01698       //std::cout <<"lij = \n" << lij << std::endl;
01699       // update pstart[j] for next time
01700       pstart[j] = startj+1;
01701       ptemp[prows[startj++]]+=lij*lij;
01702       // accumulate contributions from col j to col i in temp
01703       for (k=startj;k <degj ;k++) {
01704         rowk=prows[k];
01705         ptemp[rowk]  =ptemp[rowk]+ lij*plij[k] ;
01706         //std::cout <<"lijk = \n" << plij[k] << std::endl;
01707       }
01708     
01709       
01710       k = startj;
01711       if (k<degj) {
01712         rowk = prows[k];
01713         pcontri_cols[j] = pcontri_cols[rowk];
01714         pcontri_cols[rowk] = j;
01715        }
01716       j = nextj;
01717     }
01718     
01719     node  = &nodes[coli];
01720     prows = node->adj_list.begin();
01721     plij  = node->a.begin();
01722     diagi = pdiag[coli] - ptemp[i];
01723     //std::cout << diagi << std::endl;
01724     degi = node->adj_list.size();
01725     //         std::cout << "No of contributions to row: " << count << std::endl;
01726     //         std::cout << "No of adj elements       : " << degi   << std::endl;
01727     //         std::cout << "diagonal element: " << coli << " with value: " << pdiag[coli] << std::endl;
01728     //         std::cout << "row             : " << i    << " with value: " << ptemp[i] << std::endl;
01729     //         std::cout << "check psd:" << diagi << std::endl;
01730     //         std::cout << std::endl;
01731     if (diagi > std::max(1.e-8,1.e-10*pdiag[coli].f)) { 
01732       info_vec[0]+=log(diagi);
01733       info_vec[1]+= 1.0;
01734       //std::cout << " rank " << info_vec[1] ;
01735       diagi = sqrt(diagi);
01736       pdiag[coli] = diagi;
01737       mult = 1/diagi;
01738       // use contributions in temp to update col i
01739       for (k=0;k<degi;k++){
01740         rowk = prows[k];
01741         plij[k]-=ptemp[rowk];
01742         plij[k]*=mult;
01743         //std::cout << plij[k] << std::endl;
01744         ptemp[rowk] = 0.0;
01745       }
01746       
01747       if (degi) {
01748         rowk = prows[0];
01749         pcontri_cols[i] = pcontri_cols[rowk];
01750         pcontri_cols[rowk] = i;
01751       }
01752       
01753     }
01754     else if (diagi > -std::max(1.e-10,1.e-10*pdiag[coli].f)) {
01755       pdiag[coli] = 0;
01756       degi = node->adj_list.size();
01757       for (k=0;k<degi;k++){
01758         rowk = prows[k];
01759         plij[k] = 0.0;
01760         ptemp[rowk] = 0.0;
01761       }
01762     }
01763     else if (diagi < 0) {
01764       throw exception("SparseBGMatrix::factor matrix is not positive semi-definite");
01765     }
01766   }
01767 
01768   factor_done = 1;
01769   return 1;
01770 }

unsigned matvec::SparseBGMatrix::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 874 of file sparsebgmatrix.cpp.

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

00875 {
00876    if (path < 4) throw exception("SparseBGMatrix::factorization(path): path >= 4 is required");
00877 
00878    if (path==4) {
00879      initialize_node_list();
00880      minfilsymelm_node_list();
00881      initial_lij();
00882      factor();
00883    }
00884    if (path==6){
00885      factor();
00886    }
00887    return irank();
00888 }

BG matvec::SparseBGMatrix::getaij (  const int    ii,
const int    jj
) 
 

Return the element at (i,j) position.

Definition at line 983 of file sparsebgmatrix.cpp.

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

Referenced by initial_lij(), and operator()(). 

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

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

Definition at line 138 of file sparsebgmatrix.h.

References iiv.

Referenced by add_diag(), close(), getaij(), initialize_node_list(), insert(), q(), qTLW(), reset(), and row(). 

00138 {return iiv-1;};

int matvec::SparseBGMatrix::initial_lij (  void    ) 
 

Definition at line 1600 of file sparsebgmatrix.cpp.

References matvec::BGNodeStruct::a, matvec::BGNodeStruct::adj_list, matvec::Vector< BG >::begin(), matvec::Vector< unsigned >::begin(), matvec::Vector< BGNodeStruct >::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(). 

01601 {
01602   if (insertend) throw exception("SparseBGMatrix::initial_lij: cannot initial_lij after closing");
01603   unsigned i,j;
01604   unsigned row,col;
01605   unsigned deg,nabor;
01606   BG aij;
01607   BG *paij;
01608   unsigned *nabors;
01609   BGNodeStruct *nodes, *node;
01610 
01611   if (!elim_done) {
01612     if (!minfilsymelm_node_list()) return 0;
01613   }
01614 
01615   nodes = node_list.begin();
01616   for (i=0;i<dim;i++) {
01617     row = i+1;
01618     diag[i] = getaij(row,row);
01619     node = &nodes[i];
01620     nabors = node->adj_list.begin();
01621     deg    = node->adj_list.size();
01622 
01623     for (j=0;j<deg;j++){
01624       nabor = nabors[j];
01625       col = order(nabor);
01626       if (col < row) {
01627         aij = getaij(col,row);
01628       }
01629       else{
01630         aij = getaij(row,col);
01631       }
01632       paij = node->a.begin();
01633       paij[j] = aij;
01634     }
01635   }
01636   initial_lij_done = 1;
01637   return 1;
01638 }

int matvec::SparseBGMatrix::initialize_node_list (  void    ) 
 

Definition at line 1272 of file sparsebgmatrix.cpp.

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

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

01272                                              {
01273   int *ia   = iiv-1;
01274   int *ja   = jjv-1;
01275   BG *a = aav-1;
01276   unsigned i,j;
01277   unsigned n_rows,row,col;
01278   unsigned count;
01279 
01280   // I am going to count number of nabors for each node
01281   // I will use "mask" to do the counting!
01282 
01283   mask.resize(dim);
01284 
01285   node_list.reserve(dim);
01286   // do the counting
01287   //SDK
01288   int isrt;
01289   for(isrt=0;isrt<nonzero;isrt++) {
01290     i=srt_hash[isrt];
01291     //  for (i=1;i<=hsize;i++){
01292     //SDK
01293     if (!ia[i]) continue; //SDK
01294     row = ia[i];
01295     col = ja[i];
01296     if(row != col){
01297       mask(row) = mask(row) + 1;
01298       mask(col) = mask(col) + 1;
01299     }
01300   }
01301     // now resize the adjacent list for each node
01302   for (i=1;i<=dim;i++){
01303     node_list(i).adj_list.reserve(mask(i));
01304   }
01305   // now I am going to use mask to mark the current position in the adj_list for each node
01306   // start at position 1
01307   for (i=1;i<=dim;i++){
01308     mask(i) = 1;
01309   }
01310     // now do the filling
01311   //SDK
01312   for(isrt=0;isrt<nonzero;isrt++) {
01313     i=srt_hash[isrt];
01314     //  for (i=1;i<=hsize;i++){
01315     //SDK
01316     if (!ia[i]) continue; //SDK
01317     row = ia[i];
01318     col = ja[i];
01319     if(row!= col){
01320       node_list(row).adj_list(mask(row)) = col;
01321       node_list(col).adj_list(mask(col)) = row;
01322       mask(row) = mask(row) + 1;
01323       mask(col) = mask(col) + 1;
01324     }
01325   }
01326   initialize_node_list_done = 1;
01327   return 1;
01328 }

int matvec::SparseBGMatrix::insert (  const int    i,
const int    j,
const BG    value
) 
 

Add value to element (i,j).

Note that Indexing is from 1 instead of 0

Definition at line 933 of file sparsebgmatrix.cpp.

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

Referenced by resize(), and SparseBGMatrix(). 

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

unsigned matvec::SparseBGMatrix::irank (  void    ) 
 

Return the rank of the obhect.

Definition at line 903 of file sparsebgmatrix.cpp.

References logdet(), logdet_done, and rank.

Referenced by factorization(), and logdet(). 

00903                                {
00904   if (!logdet_done) {
00905     logdet(rank);
00906   }
00907   return rank;
00908 }

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

Definition at line 139 of file sparsebgmatrix.h.

References jjv.

Referenced by add_diag(), close(), getaij(), initialize_node_list(), insert(), q(), qTLW(), reset(), and row(). 

00139 {return jjv-1;};

BG matvec::SparseBGMatrix::logdet (  void    ) 
 

Returns log of its determinant.

Definition at line 893 of file sparsebgmatrix.cpp.

References logdet_done, logdeterm, and rank.

Referenced by irank(), and logdet(). 

00893                           {
00894   if (!logdet_done) {
00895     logdet(rank);
00896   }
00897   return logdeterm;
00898 }

BG matvec::SparseBGMatrix::logdet (  unsigned *    irank )  [inline]
 

Definition at line 131 of file sparsebgmatrix.h.

References irank(), and logdet(). 

00131 {return logdet(*irank);};

BG matvec::SparseBGMatrix::logdet (  unsigned &    irank ) 
 

Definition at line 910 of file sparsebgmatrix.cpp.

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

00910                                           {
00911    if (dim == 0) throw exception("SparseBGMatrix::logdet(): SparseBGMatrix is empty");
00912   unsigned i;
00913   BG res=0.0;
00914    if (!factor_done) {
00915      if(!factor()) return 0.0;
00916    }
00917    irank = 0;
00918    for (i=1;i<=dim;i++){
00919      if (diag(i) > 0.0000000001) {
00920        irank++;
00921        res += log(diag(i));
00922      }
00923    }
00924    logdeterm = 2*res;
00925    return (logdeterm);
00926 }

int matvec::SparseBGMatrix::minfilsymelm_node_list (  void    ) 
 

Definition at line 1330 of file sparsebgmatrix.cpp.

References a(), adjncy, matvec::Vector< BGNodeStruct >::begin(), matvec::cmp_func, matvec::comp1(), 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< BG >::reserve(), matvec::Vector< BGNodeStruct >::reserve(), matvec::Vector< unsigned >::reserve(), matvec::Vector< unsigned >::resize(), matvec::Vector< T >::resize(), matvec::Vector< BGNodeStruct >::size(), and xadj.

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

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

void matvec::SparseBGMatrix::mv (  const BG *    v,
const unsigned    n,
BG *    result
) 
 

Definition at line 439 of file sparsebgmatrix.cpp.

References aav, dim, iiv, matvec::initial_BGarray(), jjv, nonzero, and srt_hash. 

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

void matvec::SparseBGMatrix::mv (  const Vector< BG > &    v,
Vector< BG > &    result
) 
 

Definition at line 431 of file sparsebgmatrix.cpp.

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

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

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

Definition at line 122 of file sparsebgmatrix.h.

References dim. 

00122 {return dim;};

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

Definition at line 123 of file sparsebgmatrix.h.

References nonzero. 

00123 {return nonzero;};

BG matvec::SparseBGMatrix::operator() (  const int    i,
const int    j
) 
 

Retrieve element (i,j).

Definition at line 261 of file sparsebgmatrix.cpp.

References dim, and getaij(). 

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

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

Assignment operator.

Definition at line 275 of file sparsebgmatrix.cpp.

References copyfrom(). 

00276 {
00277   copyfrom(A);
00278   return *this;
00279 }

BG matvec::SparseBGMatrix::q (  const BG *    x,
const BG *    y,
const int    i1,
const int    i2
) 
 

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

SparseBGMatrix 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 595 of file sparsebgmatrix.cpp.

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

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

BG matvec::SparseBGMatrix::q (  const BG *    x,
const BG *    y
) 
 

Calculates the quadratic form x*A*y.

SparseBGMatrix 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 556 of file sparsebgmatrix.cpp.

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

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

BG matvec::SparseBGMatrix::q (  const Vector< BG > &    v1,
const Vector< BG > &    v2
) 
 

Calculates the product v1*A*v2.

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

Definition at line 542 of file sparsebgmatrix.cpp.

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

Referenced by qTLW(). 

00543 {
00544    if (dim == 0) throw exception("SparseBGMatrix is empty");
00545    if (dim != v1.size() || dim != v2.size()) throw exception("SparseBGMatrix::q(v1,v2): size not conformable");
00546    return q(v1.begin(),v2.begin());
00547 }

BG matvec::SparseBGMatrix::qTLW (  const BG *    x,
const BG *    y,
const int    i1,
const int    i2
) 
 

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

SparseBGMatrix 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 690 of file sparsebgmatrix.cpp.

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

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

BG matvec::SparseBGMatrix::qTLW (  const BG *    x,
const BG *    y
) 
 

Calculates the quadratic form x*A*y

SparseBGMatrix 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 647 of file sparsebgmatrix.cpp.

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

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

BG matvec::SparseBGMatrix::qTLW (  const Vector< BG > &    v1,
const Vector< BG > &    v2
) 
 

Calculates the quadratic form v1*A*v2.

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

Definition at line 633 of file sparsebgmatrix.cpp.

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

00634 {
00635    if (dim == 0) throw exception("SparseBGMatrix::qTLW(): it is empty");
00636    if (dim != v1.size() || dim != v2.size()) throw exception("SparseBGMatrix::qTLW(): size not conformable");
00637    return q(v1.begin(),v2.begin());
00638 }

void matvec::SparseBGMatrix::release (  void    ) 
 

Definition at line 507 of file sparsebgmatrix.cpp.

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

Referenced by ~SparseBGMatrix(). 

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

void matvec::SparseBGMatrix::release_nodestuff (  void    ) 
 

Definition at line 524 of file sparsebgmatrix.cpp.

References matvec::Vector< BG >::clear(), matvec::Vector< unsigned >::clear(), matvec::Vector< BGNodeStruct >::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(). 

00525 {
00526    node_list.clear();
00527    mask.clear();
00528    order.clear();
00529    inv_order.clear();
00530    diag.clear();
00531    initialize_node_list_done=0;
00532    elim_done=0;
00533    factor_done=0;
00534    initial_lij_done=0;
00535 }

void matvec::SparseBGMatrix::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 839 of file sparsebgmatrix.cpp.

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

00840 {
00841    if (path ==1) {
00842      initialize_node_list();
00843      minfilsymelm_node_list();
00844    }
00845    else if (path==2) {
00846      initialize_node_list();
00847      minfilsymelm_node_list();
00848      initial_lij();
00849    }
00850    else if (path==3) {
00851      initial_lij();
00852    }
00853    else if (path==4) {
00854      initialize_node_list();
00855      minfilsymelm_node_list();
00856      initial_lij();
00857   }
00858    else if (path==5) {
00859      initial_lij();
00860    }
00861 }

int matvec::SparseBGMatrix::reset (  const int    ii,
const int    jj,
const BG    value
) 
 

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

Definition at line 1024 of file sparsebgmatrix.cpp.

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

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

SparseBGMatrix & matvec::SparseBGMatrix::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 290 of file sparsebgmatrix.cpp.

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

Referenced by copyfrom(). 

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

Vector< BG > matvec::SparseBGMatrix::row (  const int    ithrow ) 
 

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

Definition at line 382 of file sparsebgmatrix.cpp.

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

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

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

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

Save the contents into a disk file.

Definition at line 1157 of file sparsebgmatrix.cpp.

References aav, dim, matvec::BG::f, iiv, and jjv. 

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

Vector< BG > matvec::SparseBGMatrix::solv (  Vector< BG >    rhs ) 
 

Definition at line 1773 of file sparsebgmatrix.cpp.

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

01774 {
01775   if (!factor_done) throw exception("SparseBGMatrix::factor must be called first");
01776   unsigned i,j;
01777   unsigned coli,rowj,corowj;
01778   unsigned degi;
01779   BG diagi,soli;
01780   Vector<BG> sol,fsol;
01781   unsigned *prows;
01782   BG *pa;
01783   BGNodeStruct *nodes, *node;
01784 
01785   sol.reserve(dim);
01786   fsol.reserve(dim);
01787 
01788   BG   *pdiag  = diag.begin();  pdiag--;
01789   BG   *psol   = sol.begin();    psol--;
01790   BG   *pfsol  = fsol.begin();   pfsol--;
01791   BG   *prhs   = rhs.begin();    prhs--;
01792   unsigned *porder = order.begin(); porder--;
01793 
01794 
01795   // forward elimination
01796   // loop for getting solution 1..dim
01797   nodes = node_list.begin(); nodes--;
01798   for (i=1;i<=dim;i++) {
01799     coli  = porder[i];
01800     diagi = pdiag[coli];
01801     node  = &nodes[coli];
01802     prows = node->adj_list.begin();
01803     pa    = node->a.begin();
01804     if (diagi == 0) { //lineraly dependent equation
01805       psol[coli] = 0.0;
01806     }
01807     else {
01808       soli = psol[coli] = prhs[coli]/diagi;
01809       //std::cout << soli << std::endl;
01810       // adjust rhs for this column
01811       // note: that diag and rhs are in original order
01812       //       but, the row subscripts of L are  in mindeg order
01813 
01814 
01815       degi = node->adj_list.size();
01816       for (j=0;j<degi;j++) {
01817         rowj = porder[prows[j]];
01818         //std::cout << pa[j] << std::endl;
01819         prhs[rowj] -= pa[j]*soli; 
01820       }
01821     }
01822   }
01823   // backward substitution
01824   // loop for getting solution dim..1
01825   for (i=dim;i>=1;i--) {
01826     coli = porder[i];
01827     diagi = pdiag[coli];
01828     node  = &nodes[coli];
01829     prows = node->adj_list.begin();
01830     pa    = node->a.begin();
01831     if (diagi == 0) { //lineraly dependent equation
01832       pfsol[coli] = 0.0;
01833     }
01834     else {
01835       // adjust this rhs for all previous solutions
01836       degi = node->adj_list.size();  
01837       soli=psol[coli];
01838       for (j=0;j<degi;j++) {
01839         rowj = porder[prows[j]];
01840         soli -=  pa[j]*pfsol[rowj];
01841       }
01842       pfsol[coli] = soli/diagi;
01843       //std::cout << pfsol[coli] << std::endl;
01844     }
01845   }
01846   //std::cout << fsol << std::endl;
01847   return fsol;
01848 }

int matvec::SparseBGMatrix::solv (  BG *    x,
const BG *    rhs
) 
 

Definition at line 1850 of file sparsebgmatrix.cpp.

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

Referenced by solve(). 

01851 {
01852   if (!factor_done) throw exception("SparseBGMatrix::factor must be called first");
01853   unsigned i,j;
01854   unsigned coli,rowj;
01855   unsigned degi;
01856   BG diagi,soli;
01857   Vector<BG> sol,temp_rhs;
01858   unsigned *prows;
01859   BG *pa;
01860   BGNodeStruct *nodes, *node;
01861 
01862   sol.reserve(dim);
01863   temp_rhs.reserve(dim);
01864   memcpy(temp_rhs.begin(),rhs,sizeof(BG)*dim);
01865 
01866   BG   *pdiag  = diag.begin();     pdiag--;
01867   BG   *psol   = sol.begin();       psol--;
01868   unsigned *porder = order.begin();    porder--;
01869   BG   *prhs   = temp_rhs.begin();  prhs--;
01870   BG   *pfsol; 
01871   pfsol = x; pfsol--;
01872  
01873 
01874   
01875   // forward elimination
01876   // loop for getting solution 1..dim
01877   nodes = node_list.begin(); nodes--;
01878   for (i=1;i<=dim;i++) {
01879     coli  = porder[i];
01880     diagi = pdiag[coli];
01881     node  = &nodes[coli];
01882     prows = node->adj_list.begin();
01883     pa    = node->a.begin();
01884 
01885     if (diagi == 0) { //lineraly dependent equation
01886       psol[coli] = 0.0;
01887     }
01888     else {
01889       soli = psol[coli] = prhs[coli]/diagi;
01890       // adjust rhs for this column
01891       // note: that diag and rhs are in original order
01892       //       but, the row subscripts of L are  in mindeg order
01893 
01894 
01895       degi = node->adj_list.size();
01896       for (j=0;j<degi;j++) {
01897         rowj = porder[prows[j]];
01898         prhs[rowj] = prhs[rowj] - pa[j]*soli;
01899       }
01900     }
01901   }
01902   // backward substitution
01903   // loop for getting solution dim..1
01904   for (i=dim;i>=1;i--) {
01905     coli = porder[i];
01906     diagi = pdiag[coli];
01907     node  = &nodes[coli];
01908     prows = node->adj_list.begin();
01909     pa    = node->a.begin();
01910     if (diagi == 0) { //lineraly dependent equation
01911       pfsol[coli] = 0.0;
01912     }
01913     else {
01914       // adjust this rhs for all previous solutions
01915       degi = node->adj_list.size();
01916       for (j=0;j<degi;j++) {
01917         rowj = porder[prows[j]];
01918         psol[coli] = psol[coli] - pa[j]*pfsol[rowj];
01919       }
01920       pfsol[coli] = psol[coli]/diagi;
01921     }
01922   }
01923   return 1;
01924 }

int matvec::SparseBGMatrix::solve (  Vector< BG > &    x,
const Vector< BG > &    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 756 of file sparsebgmatrix.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::Vector< T >::resize(), matvec::SESSION, matvec::Vector< T >::size(), and solv(). 

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

int matvec::SparseBGMatrix::solve (  BG *    x,
const BG *    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 797 of file sparsebgmatrix.cpp.

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

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

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

Definition at line 1994 of file sparsebgmatrix.cpp.

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

01994                                                                                                                                     {
01995   unsigned i,j,allzero;
01996   unsigned coli,rowj,corowj;
01997   unsigned degi;
01998   BG diagi,soli;
01999   //Vector<BG> sol,temp_rhs;
02000   unsigned *prows;
02001   BG *pa;
02002   BGNodeStruct *nodes, *node;
02003   
02004   
02005   if (!factor_done) {
02006     throw exception("SparseBGMatrix::factor must be called first");
02007     return 0;
02008   }
02009   //  sol.reserve(dim);
02010   //  temp_rhs.reserve(dim);
02011   // unsigned *row_pt=new unsigned [dim];
02012   unsigned *row_mask=row_pt;
02013   //  unsigned *col_pt=new unsigned [dim];
02014   unsigned *col_mask=col_pt;
02015   row_mask--;
02016   memset(row_pt,'\0',dim*sizeof(unsigned));
02017   col_mask--;
02018   memset(col_pt,'\0',dim*sizeof(unsigned));
02019 
02020   memcpy(temp_rhs.begin(),rhs,dim*sizeof(BG));
02021   allzero=1;
02022   BG   *pdiag  = diag.begin();     pdiag--;
02023   BG   *psol   = sol.begin();       psol--;
02024   unsigned *porder = order.begin();    porder--;
02025   unsigned *pinv_order = inv_order.begin();    pinv_order--;
02026   BG   *prhs   = temp_rhs.begin();  prhs--;
02027   BG   *pfsol; 
02028   pfsol = x; pfsol--;
02029   //memset(x,'\0',dim*sizeof(BG));
02030   initial_BGarray(x,dim,0.0);
02031 
02032   
02033   // forward elimination
02034   // loop for getting solution 1..dim
02035   nodes = node_list.begin(); nodes--;
02036   row_mask[row]=1;
02037   col_mask[row]=1;
02038   for (i=row;i<=dim;i++) {
02039     if(row_mask[i]) {
02040       coli  = porder[i];
02041       diagi = pdiag[coli];
02042       node  = &nodes[coli];
02043       prows = node->adj_list.begin();
02044       pa    = node->a.begin(); 
02045       //row_mask[i]=1;
02046       if (diagi != 0)  {
02047         soli = psol[coli] = prhs[coli]/diagi;
02048         row_mask[i]=1;
02049         // adjust rhs for this column
02050         // note: that diag and rhs are in original order
02051         //       but, the row subscripts of L are  in mindeg order
02052         
02053         
02054         degi = node->adj_list.size();
02055         for (j=0;j<degi;j++) {
02056           row_mask[prows[j]]=1;
02057           //      if(prows[j] == row) col_mask[i]=1;
02058           rowj = porder[prows[j]];
02059           prhs[rowj] -=  pa[j]*soli;
02060         }
02061       }
02062     }
02063   }
02064   // backward substitution
02065   // loop for getting solution dim..1
02066   for (i=dim;i>=row;i--) {
02067     if(row_mask[i]) {
02068       coli = porder[i];
02069       diagi = pdiag[coli];
02070       if (diagi != 0 )  {
02071         node  = &nodes[coli];
02072         prows = node->adj_list.begin();
02073         pa    = node->a.begin();
02074         // adjust this rhs for all previous solutions
02075         degi = node->adj_list.size();
02076         soli=psol[coli];
02077         for (j=0;j<degi;j++) {
02078           //row_mask[prows[j]]=1;
02079           rowj = porder[prows[j]];
02080           soli -= pa[j]*pfsol[rowj];
02081         }
02082         pfsol[coli] = soli/diagi;
02083       }
02084     }
02085   }
02086   //  delete [] row_pt;
02087   //  delete [] col_pt;
02088   return 1;
02089 }

Friends And Related Function Documentation

friend class GLMM [friend]
 

Definition at line 58 of file sparsebgmatrix.h.

Vector<BG> operator * (  const SparseBGMatrix &    A,
const Vector< BG > &    V
)  [friend]
 

SparseBGMatrix times vector operation.

Definition at line 398 of file sparsebgmatrix.cpp. 

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

SparseBGMatrix operator * (  const BG    s,
const SparseBGMatrix &    A
)  [friend]
 

Definition at line 503 of file sparsebgmatrix.cpp. 

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

SparseBGMatrix operator * (  const SparseBGMatrix &    A,
const BG    s
)  [friend]
 

Multiply SparseBGMatrix object by s.

Definition at line 479 of file sparsebgmatrix.cpp. 

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

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

Not equal operator.

Definition at line 376 of file sparsebgmatrix.cpp. 

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

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

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

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

Definition at line 355 of file sparsebgmatrix.cpp. 

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

Member Data Documentation

BG* matvec::SparseBGMatrix::aav [protected]
 

Definition at line 69 of file sparsebgmatrix.h.

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

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

Definition at line 68 of file sparsebgmatrix.h.

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

Vector<BG> matvec::SparseBGMatrix::diag [protected]
 

Definition at line 65 of file sparsebgmatrix.h.

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

unsigned matvec::SparseBGMatrix::dim [protected]
 

Definition at line 66 of file sparsebgmatrix.h.

Referenced by close(), copyfrom(), dense(), display(), display_inv_order(), display_lij(), display_node_list(), display_order(), factor(), getaij(), 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 SparseBGMatrix().

unsigned matvec::SparseBGMatrix::elim_done [protected]
 

Definition at line 73 of file sparsebgmatrix.h.

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

unsigned matvec::SparseBGMatrix::factor_done [protected]
 

Definition at line 73 of file sparsebgmatrix.h.

Referenced by copyfrom(), factor(), logdet(), release_nodestuff(), resize(), solv(), solve(), solvrow(), and SparseBGMatrix().

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

Definition at line 70 of file sparsebgmatrix.h.

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

unsigned matvec::SparseBGMatrix::hsize [protected]
 

Definition at line 66 of file sparsebgmatrix.h.

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

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

Definition at line 67 of file sparsebgmatrix.h.

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

Vector<BG> matvec::SparseBGMatrix::info_vec
 

Definition at line 80 of file sparsebgmatrix.h.

Referenced by copyfrom(), factor(), and SparseBGMatrix().

unsigned matvec::SparseBGMatrix::initial_lij_done [protected]
 

Definition at line 73 of file sparsebgmatrix.h.

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

unsigned matvec::SparseBGMatrix::initialize_node_list_done [protected]
 

Definition at line 73 of file sparsebgmatrix.h.

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

unsigned matvec::SparseBGMatrix::insertend [protected]
 

Definition at line 66 of file sparsebgmatrix.h.

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

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

Definition at line 64 of file sparsebgmatrix.h.

Referenced by display_inv_order(), minfilsymelm_node_list(), release_nodestuff(), and solvrow().

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

Definition at line 67 of file sparsebgmatrix.h.

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

Definition at line 67 of file sparsebgmatrix.h.

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

unsigned matvec::SparseBGMatrix::logdet_done [protected]
 

Definition at line 74 of file sparsebgmatrix.h.

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

BG matvec::SparseBGMatrix::logdeterm
 

Definition at line 79 of file sparsebgmatrix.h.

Referenced by logdet().

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

Definition at line 62 of file sparsebgmatrix.h.

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

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

Definition at line 62 of file sparsebgmatrix.h.

Referenced by elim_node(), and minfilsymelm_node_list().

unsigned matvec::SparseBGMatrix::max_nz [protected]
 

Definition at line 66 of file sparsebgmatrix.h.

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

Vector<BGNodeStruct> matvec::SparseBGMatrix::node_list [protected]
 

Definition at line 61 of file sparsebgmatrix.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::SparseBGMatrix::nonzero [protected]
 

Definition at line 66 of file sparsebgmatrix.h.

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

unsigned matvec::SparseBGMatrix::nsp [protected]
 

Definition at line 66 of file sparsebgmatrix.h.

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

Definition at line 63 of file sparsebgmatrix.h.

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

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

Definition at line 67 of file sparsebgmatrix.h.

unsigned matvec::SparseBGMatrix::rank [protected]
 

Definition at line 75 of file sparsebgmatrix.h.

Referenced by irank(), and logdet().

BG * matvec::SparseBGMatrix::rsp [protected]
 

Definition at line 69 of file sparsebgmatrix.h.

char matvec::SparseBGMatrix::solver_name [protected]
 

Definition at line 71 of file sparsebgmatrix.h.

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

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

Definition at line 70 of file sparsebgmatrix.h.

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

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

Definition at line 68 of file sparsebgmatrix.h.

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

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