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

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

Parent Directory Parent Directory | Revision Log Revision Log


Revision 415 - (hide annotations)
Wed Jan 4 05:37:33 2006 UTC (14 years, 1 month ago) by gross
File MIME type: text/plain
File size: 5196 byte(s)
a better way of representing the matrix format type is introduced. this is needed for the Paradiso and UMFPACK interface
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 gross 415 if (A->type & MATRIX_FORMAT_SYM) {
55 gross 411 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 gross 415 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 gross 411 } else {
66 gross 415 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 gross 411 }
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->direct==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->direct=&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->direct);
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