/[escript]/trunk/paso/src/SCSL_direct.c
ViewVC logotype

Contents of /trunk/paso/src/SCSL_direct.c

Parent Directory Parent Directory | Revision Log Revision Log


Revision 425 - (show annotations)
Tue Jan 10 04:10:39 2006 UTC (13 years, 7 months ago) by gross
File MIME type: text/plain
File size: 5196 byte(s)
The sparse solver can be called by paso now. 

the building has been change to reduce some code redundancy:
now all scons SCscripts are importing scons/esys_options.py which
imports platform specific settings. 



1 /* $Id: SCSL_direct.c 405 2005-12-22 23:05:31Z gross $ */
2
3 /**************************************************************/
4
5 /* Paso: interface to the SGI's direct solvers */
6
7 /**************************************************************/
8
9 /* Copyrights by ACcESS Australia 2003,2004,2005 */
10 /* Author: gross@access.edu.au */
11
12 /**************************************************************/
13
14 #include "Paso.h"
15 #include "SystemMatrix.h"
16 #include "SCSL.h"
17 #ifdef SCSL
18 #include <scsl_sparse.h>
19 static int TokenList[]={0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0};
20 static int Token[]={0,1,2,3,4,5,6,7,8,9,10,11,12,13,14,15,16,17,18,19};
21 static int TokenSym[]={FALSE,FALSE,FALSE,FALSE,FALSE,FALSE,FALSE,FALSE,FALSE,FALSE,FALSE,FALSE,FALSE,FALSE,FALSE,FALSE,FALSE,FALSE,FALSE,FALSE};
22 #endif
23
24 /**************************************************************/
25
26 /* free any extra stuff possibly used by the SCSL library */
27
28 void Paso_SCSL_direct_free(Paso_SystemMatrix* A) {
29 #ifdef SCSL
30 if (A->solver!=NULL) {
31 int token=*(int*)(A->solver);
32 TokenList[token]=0;
33 if (TokenSym[token]) {
34 DPSLDLT_Destroy(token);
35 } else {
36 DPSLDU_Destroy(token);
37 }
38 A->solver=NULL;
39 }
40 #endif
41 }
42 /* call the iterative solver: */
43
44 void Paso_SCSL_direct(Paso_SystemMatrix* A,
45 double* out,
46 double* in,
47 Paso_Options* options) {
48 #ifdef SCSL
49 double time0;
50 int token,l=sizeof(TokenList)/sizeof(int),reordering_method,method;
51 long long non_zeros;
52 double ops;
53
54 if (A->type & MATRIX_FORMAT_SYM) {
55 method=PASO_CHOLEVSKY;
56 } else {
57 method=PASO_DIRECT;
58 }
59 if (A->col_block_size!=1) {
60 Paso_setError(TYPE_ERROR,"Paso_SCSL_direct: linear solver can only be applied to block size 1.");
61 }
62 if (method==PASO_CHOLEVSKY) {
63 if (! (A->type & (MATRIX_FORMAT_CSC + MATRIX_FORMAT_SYM + MATRIX_FORMAT_BLK1)) )
64 Paso_setError(TYPE_ERROR,"Paso_SCSL_direct: direct solver for symmetric matrices can only be applied to symmetric CSC format with block size 1 and index offset 0.");
65 } else {
66 if (! (A->type & (MATRIX_FORMAT_CSC + MATRIX_FORMAT_BLK1)) )
67 Paso_setError(TYPE_ERROR,"Paso_SCSL_direct: direct solver can only be applied to CSC format with block size 1 and index offset 0.");
68 }
69 method=Paso_Options_getSolver(options->method,PASO_PASO,options->symmetric);
70 if (Paso_noError()) {
71
72 /* if no token has been assigned a free token must be found*/
73 if (A->solver==NULL) {
74 /* find the next available token */
75 for(token=0;token<l && TokenList[token]!=0;token++);
76 if (token==l) {
77 Paso_setError(TYPE_ERROR,"Paso_SCSL_direct: limit of number of matrices reached.");
78 } else {
79 TokenList[token] = 1;
80 TokenSym[token]=(method==PASO_CHOLEVSKY);
81 A->solver=&Token[token];
82 /* map the reordering method onto SCSL */
83 switch (options->reordering) {
84 case PASO_NO_REORDERING:
85 reordering_method=0;
86 break;
87 case PASO_MINIMUM_FILL_IN:
88 reordering_method=1;
89 break;
90 default:
91 reordering_method=2;
92 break;
93 }
94 time0=Paso_timer();
95 if (TokenSym[token]) {
96 /* DPSLDLT_Ordering(token,reordering_method); (does not work)*/
97 DPSLDLT_Preprocess(token,A->num_rows,A->pattern->ptr,A->pattern->index,&non_zeros,&ops);
98 DPSLDLT_Factor(token,A->num_rows,A->pattern->ptr,A->pattern->index,A->val);
99 if (options->verbose) printf("timing SCSL: Cholevsky factorization: %.4e sec (token = %d)\n",Paso_timer()-time0,token);
100 } else {
101 /* DPSLDU_Ordering(token,reordering_method);(does not work)*/
102 DPSLDU_Preprocess(token,A->num_rows,A->pattern->ptr,A->pattern->index,&non_zeros,&ops);
103 DPSLDU_Factor(token,A->num_rows,A->pattern->ptr,A->pattern->index,A->val);
104 if (options->verbose) printf("timing SCSL: LDU factorization: %.4e sec (token = %d)\n",Paso_timer()-time0,token);
105 }
106 }
107 }
108 }
109 time0=Paso_timer();
110
111 if (Paso_noError()) {
112 token=*(int*)(A->solver);
113 if (TokenSym[token]) {
114 DPSLDLT_Solve(token,out,in);
115 } else {
116 DPSLDU_Solve(token,out,in);
117 }
118 if (options->verbose) printf("timing SCSL: solve: %.4e sec (token = %d)\n",Paso_timer()-time0,token);
119 }
120 #else
121 Paso_setError(SYSTEM_ERROR,"Paso_SCSL_direct:SCSL is not avialble.");
122 #endif
123 }
124 /*
125 * $Log$
126 * Revision 1.2 2005/09/15 03:44:40 jgs
127 * Merge of development branch dev-02 back to main trunk on 2005-09-15
128 *
129 * Revision 1.1.2.2 2005/09/07 00:59:08 gross
130 * some inconsistent renaming fixed to make the linking work.
131 *
132 * Revision 1.1.2.1 2005/09/05 06:29:49 gross
133 * These files have been extracted from finley to define a stand alone libray for iterative
134 * linear solvers on the ALTIX. main entry through Paso_solve. this version compiles but
135 * has not been tested yet.
136 *
137 *
138 */

  ViewVC Help
Powered by ViewVC 1.1.26