/[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 1384 - (hide annotations)
Fri Jan 11 02:29:38 2008 UTC (12 years, 1 month ago) by phornby
Original Path: temp_trunk_copy/paso/src/SCSL_direct.c
File MIME type: text/plain
File size: 5875 byte(s)
Make a temp copy of the trunk before checking in the windows changes


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