## Direct sparse Cholesky solver

Das Forum zu Z88OS (Open Source) / The newsgroup for Z88OS (open source)

finzi
Import
Beiträge: 15
Registriert: Mi 28. Nov 2012, 20:10
Wohnort: Triuggio (MI) Italy

### Direct sparse Cholesky solver

Hello every body,
I'm here again with my multiple loadcases problem.
I solved the problem in teory but it is not very efficient:
- I'm testing with a brick model like this: 38000 elements, 60000 nodes (sparse matrix GS 5.5e6 Elements);
- solving with SICCG it takes 150 seconds per load case
- solving with SORCG it takes 450 seconds per load case
- solving with Cholesky is impossible for 32bit ,complied with "Pelles C" 64bit it take several hours (most probably due to 10GB memory allocation).

Then I tried to develop a Direct sparse Cholesky solver using CSparse by Timothy A. Davis.
I developped the solution in 2 steps: factorization and solution (fore - back substitution).
I had to develop also ruotines to "reload" the model, so that it's possible to get a solution without additional factorizations.

It work: factorization take 80 seconds and then 15 seconds per load case.

Here follow the code (annex. the file)
------------------------------------------------------------------------------------------------------------------
#include <cs.h>
................
.................
extern css *S ;
extern csn *N ;
/***********************************************************************
Function sp_choy_f88 - CSparse Cholesky factorization
***********************************************************************/
/***********************************************************************
This is a trial to develop a direct sparse solver using CSparse library.
Andrea Finzi <andrea.finzi@email.it>
Feb. 1, 2013
***********************************************************************/
/* CSparse Version 3.1.0 */
/* CSparse release date: Jun 1, 2012 */
/* Copyright (c) Timothy A. Davis, 2006-2012*/
int sp_choy_f88(void){
FR_INT4 i,k;
cs *T, *A ;
int order=1;
T = cs_spalloc (0, 0, 1, 1, 1) ; /* allocate result*/
if(!cs_entry(T, (csi)0, (csi)0, GS[1]))return 1; /* Fill upper trinagular part of spars matrix T*/
for(i=2;i<=nfg;i++){ /* r=1 made manually because iez[ip[r-1]+1] is not accesibe*/
for(k=ip[i-1]+1;k<=ip;k++)/* T is zero-based and A is one-based*/
if (!cs_entry(T, (csi)(iez[k]-1), (csi)(i-1), GS[k]))return 1;
/* T(c-1,r-1)=A(r,c) (r=1,nfg; c=iez[ip[r-1]+1],iez[ip[r]]) */
}
A = cs_compress (T) ; /* A = compressed-column form of T */
cs_spfree (T) ; /* clear T */
S = cs_schol (order, A) ; /* ordering and symbolic analysis */
N = cs_chol (A, S) ; /* numeric Cholesky factorization */
cs_spfree (A) ; /* clear A */
return 0;
}

/***********************************************************************
Function sp_choy_s88 - CSparse Cholesky solution
***********************************************************************/
/***********************************************************************
This is a trial to develop a direct sparse solver.
Andrea Finzi <andrea.finzi@email.it>
Feb. 1, 2013
***********************************************************************/
/* CSparse Version 3.1.0*/
/* CSparse release date: Jun 1, 2012*/
/* Copyright (c) Timothy A. Davis, 2006-2012*/
int sp_choy_s88(void){
FR_INT4 i;
double *x, *b ;

x = cs_malloc (nfg, sizeof (double)) ; /* get workspace */
b = cs_malloc (nfg, sizeof (double)) ; /* get workspace */
for(i= 1;i <= nfg;i++) b[i-1] = rs;
if (S && N && x && b) {
cs_ipvec (S->pinv, b, x, nfg) ; /* x = P*b */
cs_lsolve (N->L, x) ; /* x = L\x */
cs_ltsolve (N->L, x) ; /* x = L'\x */
cs_pvec (S->pinv, x, b, nfg) ; /* b = P'*x */
}
for(i= 1;i <= nfg;i++) rs = b[i-1];
cs_free (x) ;
cs_free (b) ;
return 0;
}
----------------------------------------------------------------------------------------------------------
Dateianhänge
choy88s.zip
Direct sparse Cholesky solver
Think positive, if you can!

auroraIco