/[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 1811 - (show annotations)
Thu Sep 25 23:11:13 2008 UTC (10 years, 8 months ago) by ksteube
File MIME type: text/plain
File size: 5840 byte(s)
Copyright updated in all files

1
2 /*******************************************************
3 *
4 * Copyright (c) 2003-2008 by University of Queensland
5 * Earth Systems Science Computational Center (ESSCC)
6 * http://www.uq.edu.au/esscc
7 *
8 * Primary Business: Queensland, Australia
9 * Licensed under the Open Software License version 3.0
10 * http://www.opensource.org/licenses/osl-3.0.php
11 *
12 *******************************************************/
13
14
15 /**************************************************************/
16
17 /* Paso: interface to the SGI's direct solvers */
18
19 /**************************************************************/
20
21 /* Copyrights by ACcESS Australia 2003,2004,2005 */
22 /* Author: gross@access.edu.au */
23
24 /**************************************************************/
25
26 #include "Paso.h"
27 #include "SystemMatrix.h"
28 #include "SCSL.h"
29 #ifdef SCSL
30 #include <scsl_sparse.h>
31 static int TokenList[]={0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0};
32 static int Token[]={0,1,2,3,4,5,6,7,8,9,10,11,12,13,14,15,16,17,18,19};
33 static int TokenSym[]={FALSE,FALSE,FALSE,FALSE,FALSE,FALSE,FALSE,FALSE,FALSE,FALSE,FALSE,FALSE,FALSE,FALSE,FALSE,FALSE,FALSE,FALSE,FALSE,FALSE};
34 #endif
35
36 /**************************************************************/
37
38 /* free any extra stuff possibly used by the SCSL library */
39
40 void Paso_SCSL_direct_free(Paso_SystemMatrix* A) {
41 #ifdef SCSL
42 if (A->solver!=NULL) {
43 int token=*(int*)(A->solver);
44 TokenList[token]=0;
45 if (TokenSym[token]) {
46 DPSLDLT_Destroy(token);
47 } else {
48 DPSLDU_Destroy(token);
49 }
50 A->solver=NULL;
51 }
52 #endif
53 }
54 /* call the iterative solver: */
55
56 void Paso_SCSL_direct(Paso_SystemMatrix* A,
57 double* out,
58 double* in,
59 Paso_Options* options,
60 Paso_Performance* pp) {
61 #ifdef SCSL
62 double time0;
63 int token,l=sizeof(TokenList)/sizeof(int),reordering_method,method;
64 long long non_zeros;
65 double ops;
66
67 if (A->type & MATRIX_FORMAT_SYM) {
68 method=PASO_CHOLEVSKY;
69 } else {
70 method=PASO_DIRECT;
71 }
72 if (A->col_block_size!=1) {
73 Paso_setError(TYPE_ERROR,"Paso_SCSL_direct: linear solver can only be applied to block size 1.");
74 }
75 if (method==PASO_CHOLEVSKY) {
76 if (! (A->type & (MATRIX_FORMAT_CSC + MATRIX_FORMAT_SYM + MATRIX_FORMAT_BLK1)) )
77 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.");
78 } else {
79 if (! (A->type & (MATRIX_FORMAT_CSC + MATRIX_FORMAT_BLK1)) )
80 Paso_setError(TYPE_ERROR,"Paso_SCSL_direct: direct solver can only be applied to CSC format with block size 1 and index offset 0.");
81 }
82 method=Paso_Options_getSolver(options->method,PASO_PASO,options->symmetric);
83 Performance_startMonitor(pp,PERFORMANCE_ALL);
84 if (Paso_noError()) {
85
86 /* if no token has been assigned a free token must be found*/
87 if (A->solver==NULL) {
88 /* find the next available token */
89 for(token=0;token<l && TokenList[token]!=0;token++);
90 if (token==l) {
91 Paso_setError(TYPE_ERROR,"Paso_SCSL_direct: limit of number of matrices reached.");
92 } else {
93 TokenList[token] = 1;
94 TokenSym[token]=(method==PASO_CHOLEVSKY);
95 A->solver=&Token[token];
96 /* map the reordering method onto SCSL */
97 switch (options->reordering) {
98 case PASO_NO_REORDERING:
99 reordering_method=0;
100 break;
101 case PASO_MINIMUM_FILL_IN:
102 reordering_method=1;
103 break;
104 default:
105 reordering_method=2;
106 break;
107 }
108 time0=Paso_timer();
109 if (TokenSym[token]) {
110 /* DPSLDLT_Ordering(token,reordering_method); (does not work)*/
111 DPSLDLT_Preprocess(token,A->mainBlock->numRows,A->mainBlock->pattern->ptr,A->mainBlock->pattern->index,&non_zeros,&ops);
112 DPSLDLT_Factor(token,A->mainBlock->numRows,A->mainBlock->pattern->ptr,A->mainBlock->pattern->index,A->mainBlock->val);
113 if (options->verbose) printf("timing SCSL: Cholevsky factorization: %.4e sec (token = %d)\n",Paso_timer()-time0,token);
114 } else {
115 /* DPSLDU_Ordering(token,reordering_method);(does not work)*/
116 DPSLDU_Preprocess(token,A->mainBlock->numRows,A->mainBlock->pattern->ptr,A->mainBlock->pattern->index,&non_zeros,&ops);
117 DPSLDU_Factor(token,A->mainBlock->numRows,A->mainBlock->pattern->ptr,A->mainBlock->pattern->index,A->mainBlock->val);
118 if (options->verbose) printf("timing SCSL: LDU factorization: %.4e sec (token = %d)\n",Paso_timer()-time0,token);
119 }
120 }
121 }
122 }
123 time0=Paso_timer();
124 if (Paso_noError()) {
125 token=*(int*)(A->solver);
126 if (TokenSym[token]) {
127 DPSLDLT_Solve(token,out,in);
128 } else {
129 DPSLDU_Solve(token,out,in);
130 }
131 if (options->verbose) printf("timing SCSL: solve: %.4e sec (token = %d)\n",Paso_timer()-time0,token);
132 }
133 Performance_stopMonitor(pp,PERFORMANCE_ALL);
134 #else
135 Paso_setError(SYSTEM_ERROR,"Paso_SCSL_direct:SCSL is not avialble.");
136 #endif
137 }
138 /*
139 * $Log$
140 * Revision 1.2 2005/09/15 03:44:40 jgs
141 * Merge of development branch dev-02 back to main trunk on 2005-09-15
142 *
143 * Revision 1.1.2.2 2005/09/07 00:59:08 gross
144 * some inconsistent renaming fixed to make the linking work.
145 *
146 * Revision 1.1.2.1 2005/09/05 06:29:49 gross
147 * These files have been extracted from finley to define a stand alone libray for iterative
148 * linear solvers on the ALTIX. main entry through Paso_solve. this version compiles but
149 * has not been tested yet.
150 *
151 *
152 */

  ViewVC Help
Powered by ViewVC 1.1.26