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

  ViewVC Help
Powered by ViewVC 1.1.26