472,791 Members | 1,878 Online

# A fast column based Cholesky factorization in C++ needed

Hi all,

I am currently coding a sparse factorization program for my project. I
intend to make a comparison of the column based method and multifrontal
method in terms of running time, but I am not sure if my implementation of
the column based factorization is doing the method just. So I am looking for
some other source of implementation in C++ to see if there is any possible
optimization that I have missed out. Anyone know of such source code
available? Any reference and link is welcome and would be greatly
appreciated. Thanks.

Thuan Seah
Jul 22 '05 #1
4 2909
Tan Thuan Seah wrote:
I am currently coding a sparse factorization program for my project. I
intend to make a comparison of the column based method and multifrontal
method in terms of running time, but I am not sure if my implementation of
the column based factorization is doing the method just. So I am looking for
some other source of implementation in C++ to see if there is any possible
optimization that I have missed out.
Anyone know of such source code
available?
Try comp.sources.wanted.
Any reference and link is welcome and would be greatly
appreciated. Thanks.

Jul 22 '05 #2
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;
Real Ljk=0.0;
Real Lik=0.0;
Real Lij=0.0;
Real *update;
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];
vector<int> blank;
blank.clear();
for(int i=0; i<ncol; i++) {
update[i] = 0.0f;
}
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
int k = *p1;
Ljk = L[k][counter];
counter++;
// update the diagonal
update[j] += Ljk*Ljk;
// cmod(j,k) go through all the off diagonal elements of
Lk
// update
update[*p2] += L[k][counter]*Ljk;
counter++;
}
}
/********** 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)
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!
Jul 22 '05 #3
"Tan Thuan Seah" <u2******@anu.edu.au> wrote in message news:<41********@clarion.carno.net.au>...
Hi all,

I am currently coding a sparse factorization program for my project. I
intend to make a comparison of the column based method and multifrontal
method in terms of running time, but I am not sure if my implementation of
the column based factorization is doing the method just. So I am looking for
some other source of implementation in C++ to see if there is any possible
optimization that I have missed out. Anyone know of such source code
available? Any reference and link is welcome and would be greatly
appreciated. Thanks.

try
http://gams.nist.gov/
and
http://math.nist.gov/
Jul 22 '05 #4
Tan Thuan Seah wrote:
I am currently coding a sparse factorization program for my project.
I intend to make a comparison of the column based method
and multifrontal method in terms of running time,
but I am not sure if my implementation of the column based factorization
is doing the method just.
So I am looking for some other source of implementation in C++
to see if there is any possible optimization that I have missed out.
Anyone know of such source code available?
Any reference and link is welcome and would be greatly appreciated.
Take a look at
The C++ Scalar, Vector, Matrix and Tensor class library

http://www.netwood.net/~edwin/svmtl/
cat svmtl/src/square/Square.ccP

Jul 22 '05 #5

This thread has been closed and replies have been disabled. Please start a new discussion.