/[escript]/trunk/esys2/finley/src/finleyC/SCSL/SCSL_direct.c
ViewVC logotype

Contents of /trunk/esys2/finley/src/finleyC/SCSL/SCSL_direct.c

Parent Directory Parent Directory | Revision Log Revision Log


Revision 126 - (show annotations)
Fri Jul 22 03:53:08 2005 UTC (13 years, 8 months ago) by jgs
File MIME type: text/plain
File size: 5075 byte(s)
Merge of development branch back to main trunk on 2005-07-22

1 /* $Id$ */
2
3 /**************************************************************/
4
5 /* Finley: interface to the direct solvers */
6
7 /**************************************************************/
8
9 /* Copyrights by ACcESS Australia 2003 */
10 /* Author: gross@access.edu.au */
11
12 /**************************************************************/
13
14 #include "Finley.h"
15 #include "System.h"
16 #include "SCSL.h"
17 #if DIRECT_SOLVER == SGI_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 Finley_SCSL_direct_free(Finley_SystemMatrix* A) {
29 #if DIRECT_SOLVER == SGI_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 Finley_SCSL_direct(Finley_SystemMatrix* A,
45 double* out,
46 double* in,
47 Finley_SolverOptions* options) {
48 #if DIRECT_SOLVER == SGI_SCSL
49 double time0;
50 int token,l=sizeof(TokenList)/sizeof(int),reordering_method,method;
51 long long non_zeros;
52 double ops;
53 if (options->symmetric) {
54 method=ESCRIPT_CHOLEVSKY;
55 } else {
56 method=ESCRIPT_DIRECT;
57 }
58 if (A->col_block_size!=1) {
59 Finley_ErrorCode=TYPE_ERROR;
60 sprintf(Finley_ErrorMsg,"SCSL linear solver can only be applied to block size 1.");
61 return;
62 }
63 if (method==ESCRIPT_CHOLEVSKY) {
64 if (A->type!=CSC_SYM) {
65 Finley_ErrorCode=TYPE_ERROR;
66 sprintf(Finley_ErrorMsg,"SGI SCSL direct solver for symmetric matrices can only be applied to symmetric CSC format.");
67 return;
68 }
69 } else {
70 if (A->type!=CSC) {
71 Finley_ErrorCode=TYPE_ERROR;
72 sprintf(Finley_ErrorMsg,"SGI SCSL direct solver can only be applied to CSC format.");
73 return;
74 }
75 }
76 /* if no token has been assigned a free token must be found*/
77 if (A->direct==NULL) {
78 /* find the next available token */
79 for(token=0;token<l && TokenList[token]!=0;token++);
80 if (token==l) {
81 Finley_ErrorCode=TYPE_ERROR;
82 sprintf(Finley_ErrorMsg,"SCSL linear solver can only deal with %d matrices simultaneously.",l);
83 return;
84 } else {
85 TokenList[token] = 1;
86 TokenSym[token]=(method==ESCRIPT_CHOLEVSKY);
87 A->direct=&Token[token];
88 /* map the reordering method onto SCSL */
89 switch (options->reordering) {
90 case ESCRIPT_NO_REORDERING:
91 reordering_method=0;
92 break;
93 case ESCRIPT_MINIMUM_FILL_IN:
94 reordering_method=1;
95 break;
96 default:
97 reordering_method=2;
98 break;
99 }
100 time0=Finley_timer();
101 if (TokenSym[token]) {
102 /* DPSLDLT_Ordering(token,reordering_method); (does not work)*/
103 DPSLDLT_Preprocess(token,A->num_rows,A->pattern->ptr,A->pattern->index,&non_zeros,&ops);
104 DPSLDLT_Factor(token,A->num_rows,A->pattern->ptr,A->pattern->index,A->val);
105 if (options->verbose) printf("timing SCSL: Cholevsky factorization: %.4e sec (token = %d)\n",Finley_timer()-time0,token);
106 } else {
107 /* DPSLDU_Ordering(token,reordering_method);(does not work)*/
108 DPSLDU_Preprocess(token,A->num_rows,A->pattern->ptr,A->pattern->index,&non_zeros,&ops);
109 DPSLDU_Factor(token,A->num_rows,A->pattern->ptr,A->pattern->index,A->val);
110 if (options->verbose) printf("timing SCSL: LDU factorization: %.4e sec (token = %d)\n",Finley_timer()-time0,token);
111 }
112 }
113 }
114 time0=Finley_timer();
115 token=*(int*)(A->direct);
116 if (TokenSym[token]) {
117 DPSLDLT_Solve(token,out,in);
118 } else {
119 DPSLDU_Solve(token,out,in);
120 }
121 if (options->verbose) printf("timing SCSL: solve: %.4e sec (token = %d)\n",Finley_timer()-time0,token);
122 #else
123 Finley_ErrorCode=SYSTEM_ERROR;
124 sprintf(Finley_ErrorMsg,"No SCSL direct solver available.");
125 #endif
126 }
127 /*
128 * $Log$
129 * Revision 1.3 2005/07/22 03:53:08 jgs
130 * Merge of development branch back to main trunk on 2005-07-22
131 *
132 * Revision 1.2 2004/12/15 07:08:34 jgs
133 * *** empty log message ***
134 *
135 * Revision 1.1.2.2 2005/07/18 10:34:55 gross
136 * some informance improvements when reading meshes
137 *
138 * Revision 1.1.2.1 2004/11/12 06:58:20 gross
139 * a lot of changes to get the linearPDE class running: most important change is that there is no matrix format exposed to the user anymore. the format is chosen by the Domain according to the solver and symmetry
140 *
141 *
142 */

Properties

Name Value
svn:eol-style native
svn:keywords Author Date Id Revision

  ViewVC Help
Powered by ViewVC 1.1.26