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

Annotation of /trunk-mpi-branch/paso/src/SCSL_direct.c

Parent Directory Parent Directory | Revision Log Revision Log


Revision 411 - (hide annotations)
Tue Jan 3 00:23:48 2006 UTC (13 years, 10 months ago) by gross
Original Path: trunk/paso/src/SCSL_direct.c
File MIME type: text/plain
File size: 5029 byte(s)
SCSL interface has moved

1 gross 411 /* $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->direct!=NULL) {
31     int token=*(int*)(A->direct);
32     TokenList[token]=0;
33     if (TokenSym[token]) {
34     DPSLDLT_Destroy(token);
35     } else {
36     DPSLDU_Destroy(token);
37     }
38     A->direct=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==CSC_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!=CSC_SYM) {
64     Paso_setError(TYPE_ERROR,"Paso_SCSL_direct: direct solver for symmetric matrices can only be applied to symmetric CSC format.");
65     }
66     } else {
67     if (A->type!=CSC) {
68     Paso_setError(TYPE_ERROR,"Paso_SCSL_direct: direct solver can only be applied to CSC format.");
69     }
70     }
71     method=Paso_Options_getSolver(options->method,PASO_PASO,options->symmetric);
72     if (Paso_noError()) {
73    
74     /* if no token has been assigned a free token must be found*/
75     if (A->direct==NULL) {
76     /* find the next available token */
77     for(token=0;token<l && TokenList[token]!=0;token++);
78     if (token==l) {
79     Paso_setError(TYPE_ERROR,"Paso_SCSL_direct: limit of number of matrices reached.");
80     } else {
81     TokenList[token] = 1;
82     TokenSym[token]=(method==PASO_CHOLEVSKY);
83     A->direct=&Token[token];
84     /* map the reordering method onto SCSL */
85     switch (options->reordering) {
86     case PASO_NO_REORDERING:
87     reordering_method=0;
88     break;
89     case PASO_MINIMUM_FILL_IN:
90     reordering_method=1;
91     break;
92     default:
93     reordering_method=2;
94     break;
95     }
96     time0=Paso_timer();
97     if (TokenSym[token]) {
98     /* DPSLDLT_Ordering(token,reordering_method); (does not work)*/
99     DPSLDLT_Preprocess(token,A->num_rows,A->pattern->ptr,A->pattern->index,&non_zeros,&ops);
100     DPSLDLT_Factor(token,A->num_rows,A->pattern->ptr,A->pattern->index,A->val);
101     if (options->verbose) printf("timing SCSL: Cholevsky factorization: %.4e sec (token = %d)\n",Paso_timer()-time0,token);
102     } else {
103     /* DPSLDU_Ordering(token,reordering_method);(does not work)*/
104     DPSLDU_Preprocess(token,A->num_rows,A->pattern->ptr,A->pattern->index,&non_zeros,&ops);
105     DPSLDU_Factor(token,A->num_rows,A->pattern->ptr,A->pattern->index,A->val);
106     if (options->verbose) printf("timing SCSL: LDU factorization: %.4e sec (token = %d)\n",Paso_timer()-time0,token);
107     }
108     }
109     }
110     }
111     time0=Paso_timer();
112    
113     if (Paso_noError()) {
114     token=*(int*)(A->direct);
115     if (TokenSym[token]) {
116     DPSLDLT_Solve(token,out,in);
117     } else {
118     DPSLDU_Solve(token,out,in);
119     }
120     if (options->verbose) printf("timing SCSL: solve: %.4e sec (token = %d)\n",Paso_timer()-time0,token);
121     }
122     #else
123     Paso_setError(SYSTEM_ERROR,"Paso_SCSL_direct:SCSL is not avialble.");
124     #endif
125     }
126     /*
127     * $Log$
128     * Revision 1.2 2005/09/15 03:44:40 jgs
129     * Merge of development branch dev-02 back to main trunk on 2005-09-15
130     *
131     * Revision 1.1.2.2 2005/09/07 00:59:08 gross
132     * some inconsistent renaming fixed to make the linking work.
133     *
134     * Revision 1.1.2.1 2005/09/05 06:29:49 gross
135     * These files have been extracted from finley to define a stand alone libray for iterative
136     * linear solvers on the ALTIX. main entry through Paso_solve. this version compiles but
137     * has not been tested yet.
138     *
139     *
140     */

  ViewVC Help
Powered by ViewVC 1.1.26