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 ----