/[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 584 - (hide annotations)
Thu Mar 9 23:03:38 2006 UTC (13 years, 11 months ago) by gross
File MIME type: text/plain
File size: 5344 byte(s)
eigenvalues: compiles and passes tests on altix now
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 gross 425 if (A->solver!=NULL) {
31     int token=*(int*)(A->solver);
32 gross 411 TokenList[token]=0;
33     if (TokenSym[token]) {
34     DPSLDLT_Destroy(token);
35     } else {
36     DPSLDU_Destroy(token);
37     }
38 gross 425 A->solver=NULL;
39 gross 411 }
40     #endif
41     }
42     /* call the iterative solver: */
43    
44     void Paso_SCSL_direct(Paso_SystemMatrix* A,
45     double* out,
46     double* in,
47 gross 584 Paso_Options* options,
48     Paso_Performance* pp) {
49 gross 411 #ifdef SCSL
50     double time0;
51     int token,l=sizeof(TokenList)/sizeof(int),reordering_method,method;
52     long long non_zeros;
53     double ops;
54    
55 gross 415 if (A->type & MATRIX_FORMAT_SYM) {
56 gross 411 method=PASO_CHOLEVSKY;
57     } else {
58     method=PASO_DIRECT;
59     }
60     if (A->col_block_size!=1) {
61     Paso_setError(TYPE_ERROR,"Paso_SCSL_direct: linear solver can only be applied to block size 1.");
62     }
63     if (method==PASO_CHOLEVSKY) {
64 gross 415 if (! (A->type & (MATRIX_FORMAT_CSC + MATRIX_FORMAT_SYM + MATRIX_FORMAT_BLK1)) )
65     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.");
66 gross 411 } else {
67 gross 415 if (! (A->type & (MATRIX_FORMAT_CSC + MATRIX_FORMAT_BLK1)) )
68     Paso_setError(TYPE_ERROR,"Paso_SCSL_direct: direct solver can only be applied to CSC format with block size 1 and index offset 0.");
69 gross 411 }
70     method=Paso_Options_getSolver(options->method,PASO_PASO,options->symmetric);
71 gross 584 Performance_startMonitor(pp,PERFORMANCE_ALL);
72 gross 411 if (Paso_noError()) {
73    
74     /* if no token has been assigned a free token must be found*/
75 gross 425 if (A->solver==NULL) {
76 gross 411 /* 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 gross 425 A->solver=&Token[token];
84 gross 411 /* 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     if (Paso_noError()) {
113 gross 425 token=*(int*)(A->solver);
114 gross 411 if (TokenSym[token]) {
115     DPSLDLT_Solve(token,out,in);
116     } else {
117     DPSLDU_Solve(token,out,in);
118     }
119     if (options->verbose) printf("timing SCSL: solve: %.4e sec (token = %d)\n",Paso_timer()-time0,token);
120     }
121 gross 584 Performance_stopMonitor(pp,PERFORMANCE_ALL);
122 gross 411 #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