Below is my code. I am not sure how readable is it, hope someone can suggest
something. Thanks.
Thuan Seah
INPUT:
Li The non zero indices for the Cholesky factor
L The storage for the L. On input, it hold the non zero
entries of matrix A and fill in entries has value of 0.0.
col The number of column matrix A has.
void ColFact(vector<vector<int> > &Li, int ncol, vector<vector<Real> > &L) {
/*******Initialization***************************** ********/
// machine epsilon
const double epsilon = FLT_MIN*100.0;
int index=0;
int counter=0;
int *head;
Real Ljk=0.0;
Real Lik=0.0;
Real Lij=0.0;
Real *update;
vector<vector<int> > link;
vector<int>::iterator indexstart;
vector<int>::iterator offset;
vector<int>::iterator valstart;
vector<int>::iterator p1;
vector<int>::iterator p2;
vector<Real>::iterator dp1;
update = new Real[ncol];
head = new int[ncol];
link.reserve(ncol);
vector<int> blank;
blank.clear();
for(int i=0; i<ncol; i++) {
update[i] = 0.0f;
head[i] = 1;
link.push_back(blank);
}
cout << "STAGE 3: NUMERIC FACTORIZATION USING COLUMN SPARSE
FACTORIZATION\n";
/**** End of Initialization************************************ ******/
/****Start of Column factorization***********************************/
// outer loop
for(int j=0; j<ncol; j++) {
// inner loop
for(p1=link[j].begin(); p1!=link[j].end(); ++p1) {
int k = *p1;
counter=head[k];
Ljk = L[k][counter];
counter++;
// update the diagonal
update[j] += Ljk*Ljk;
// cmod(j,k) go through all the off diagonal elements of
Lk
for(p2=Li[k].begin()+head[k]; p2!=Li[k].end(); ++p2) {
// update
update[*p2] += L[k][counter]*Ljk;
counter++;
}
head[k] = head[k]+1;
}
/********** Error detection ***************************/
L[j][0] = L[j][0] - update[j];
if(abs(L[j][0]) < epsilon) {
cout << "ERROR: ONE OF THE DIAGONAL ENTRIES IN L IS ZERO\n";
cout << "PROGRAM EXITING..." << endl;
}
if(L[j][0] < 0.0f) {
cout << "ERROR: DIAGONAL ENTRIES IN L" << j << " IS
NEGATIVE\n";
cout << "PROGRAM EXITING..." << endl;
exit(1);
}
/******* Finalize computation for the column **************/
L[j][0] = sqrt(L[j][0]);
update[j] = 0.0f;
counter=1;
for(p1=Li[j].begin(); p1!=Li[j].end(); ++p1) {
L[j][counter] = (L[j][counter]-update[*p1])/L[j][0];
if(abs(L[j][counter]) > epsilon)
link[*p1].push_back(j);
update[*p1] = 0.0f;
counter++;
}
}
cout << "STAGE 3: NUMERIC FACTORIZATION DONE\n";
/******* End of Column Factorization**********************************/
}
------------ And now a word from our sponsor ------------------
Do your users want the best web-email gateway? Don't let your
customers drift off to free webmail services install your own
web gateway!
-- See
http://netwinsite.com/sponsor/sponsor_webmail.htm ----
