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

1 gross 411
2 ksteube 1312 /*******************************************************
3 ksteube 1811 *
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 dhawcroft 631
14 ksteube 1811
15 gross 411 /**************************************************************/
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 gross 425 if (A->solver!=NULL) {
43     int token=*(int*)(A->solver);
44 gross 411 TokenList[token]=0;
45     if (TokenSym[token]) {
46     DPSLDLT_Destroy(token);
47     } else {
48     DPSLDU_Destroy(token);
49     }
50 gross 425 A->solver=NULL;
51 gross 411 }
52     #endif
53     }
54     /* call the iterative solver: */
55    
56     void Paso_SCSL_direct(Paso_SystemMatrix* A,
57     double* out,
58     double* in,
59 gross 584 Paso_Options* options,
60     Paso_Performance* pp) {
61 gross 411 #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 gross 415 if (A->type & MATRIX_FORMAT_SYM) {
68 gross 411 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 gross 415 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 gross 411 } else {
79 gross 415 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 gross 411 }
82     method=Paso_Options_getSolver(options->method,PASO_PASO,options->symmetric);
83 gross 584 Performance_startMonitor(pp,PERFORMANCE_ALL);
84 gross 411 if (Paso_noError()) {
85    
86     /* if no token has been assigned a free token must be found*/
87 gross 425 if (A->solver==NULL) {
88 gross 411 /* 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 gross 425 A->solver=&Token[token];
96 gross 411 /* 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 ksteube 1312 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 gross 411 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 ksteube 1312 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 gross 411 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 gross 425 token=*(int*)(A->solver);
126 gross 411 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 gross 584 Performance_stopMonitor(pp,PERFORMANCE_ALL);
134 gross 411 #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