/[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 411 - (show annotations)
Tue Jan 3 00:23:48 2006 UTC (13 years, 6 months ago) by gross
File MIME type: text/plain
File size: 5029 byte(s)
SCSL interface has moved

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->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