/[escript]/trunk/paso/src/SystemMatrix_saveHB.c
ViewVC logotype

Annotation of /trunk/paso/src/SystemMatrix_saveHB.c

Parent Directory Parent Directory | Revision Log Revision Log


Revision 686 - (hide annotations)
Mon Mar 27 22:17:16 2006 UTC (16 years, 8 months ago) by gross
File MIME type: text/plain
File size: 7267 byte(s)
small fixes
1 jgs 150 /* $Id$ */
2    
3 dhawcroft 631 /*
4     ********************************************************************************
5 dhawcroft 633 * Copyright 2006 by ACcESS MNRF *
6 dhawcroft 631 * *
7     * http://www.access.edu.au *
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 jgs 150 /**************************************************************/
15    
16     /* Paso: SystemMatrix is saved to Harwell-Boeing format */
17    
18     /**************************************************************/
19    
20     /* Copyright: ACcESS Australia 2005 */
21     /* Author: imran@esscc.uq.edu.au */
22    
23     /**************************************************************/
24    
25     #include "Paso.h"
26     #include "SystemMatrix.h"
27    
28     /* TODO: Refactor the stuff in here into hbio, just like mmio! */
29    
30     static dim_t M, N, nz;
31    
32     static int calc_digits( int );
33     static void fmt_str( int, int, int*, int*, int*, char*, char* );
34     static void print_data( FILE*, int, int, int, char*, void*, int, int );
35     static void generate_HB( FILE*, dim_t*, dim_t*, double* );
36    
37     /* function to get number of digits in an integer */
38     int calc_digits( int var )
39     {
40     int digits = 1;
41     while( (var/=10) )
42     digits++;
43    
44     return digits;
45     }
46    
47     /* function to generate the format string.
48     *
49     * use maxlen to determine no. of entries per line
50     * use nvalues to determine no. of lines
51     */
52     void fmt_str( int nvalues, int integer, int *width, int *nlines, int *nperline, char *pfmt, char *fmt )
53     {
54     int per_line;
55     int maxlen = *width;
56    
57     if( integer && maxlen < 10 )
58     maxlen = 10;
59     else
60     maxlen = 13;
61    
62     per_line = 80 / maxlen;
63     *nlines = nvalues / per_line;
64     if( nvalues % per_line )
65     (*nlines)++;
66    
67     *nperline = per_line;
68     if( integer )
69     {
70     sprintf( pfmt, "(%dI%d)", per_line, maxlen );
71     sprintf( fmt, "%%%dd", maxlen );
72     }
73     else
74     {
75     sprintf( pfmt, "(1P%dE%d.6)", per_line, maxlen );
76     sprintf( fmt, "%%%d.6E", maxlen );
77     }
78     *width = maxlen;
79     }
80    
81     /* function to print the actual data in the right format */
82     void print_data( FILE *fp, int n_perline, int width, int nval, char *fmt, void *ptr, int integer, int adjust )
83     {
84 gross 686 double *data = ptr;
85 jgs 150 int entries_done = 0;
86 gross 686 int padding, i, j;
87 jgs 150 char pad_fmt[10];
88    
89     padding = 80 - n_perline*width;
90     sprintf( pad_fmt, "%%%dc", padding );
91    
92     if( adjust != 1 )
93     adjust = 0;
94    
95     if( integer )
96     {
97     dim_t *data = ptr;
98 gross 686 for(i=0; i<nval; i++ )
99 jgs 150 {
100     fprintf( fp, fmt, data[i]+adjust );
101     entries_done++;
102     if( entries_done == n_perline )
103     {
104     if( padding )
105     fprintf( fp, pad_fmt, ' ' );
106     fprintf( fp, "\n" );
107     entries_done = 0;
108     }
109     }
110     }
111     else
112     {
113 gross 686 for(i=0; i<nval; i++ )
114 jgs 150 {
115     fprintf( fp, fmt, data[i] );
116     entries_done++;
117     if( entries_done == n_perline )
118     {
119     if( padding )
120     fprintf( fp, pad_fmt, ' ' );
121     fprintf( fp, "\n" );
122     entries_done = 0;
123     }
124     }
125     }
126     if( entries_done )
127     {
128     sprintf( pad_fmt, "%%%dc\n", (80 - entries_done*width) );
129     fprintf( fp, pad_fmt, ' ' );
130     }
131     }
132    
133     void generate_HB( FILE *fp, dim_t *col_ptr, dim_t *row_ind, double *val )
134     {
135     char buffer[81];
136    
137     int val_lines, ind_lines, ptr_lines;
138     int val_perline, ind_perline, ptr_perline;
139     int val_width, ind_width, ptr_width;
140     char ptr_pfmt[7], ind_pfmt[7], val_pfmt[11];
141     char ptr_fmt[10], ind_fmt[10], val_fmt[10];
142    
143     /* line 1 */
144     sprintf( buffer, "%-72s%-8s", "Matrix Title", "Key" );
145     buffer[80] = '\0';
146     fprintf( fp, "%s\n", buffer );
147    
148     /* line 2 */
149     ptr_width = calc_digits( nz+1 );
150     fmt_str( N+1, 1, &ptr_width, &ptr_lines, &ptr_perline, ptr_pfmt, ptr_fmt );
151     ind_width = calc_digits( N );
152     fmt_str( nz, 1, &ind_width, &ind_lines, &ind_perline, ind_pfmt, ind_fmt );
153     val_width = 13;
154     fmt_str( nz, 0, &val_width, &val_lines, &val_perline, val_pfmt, val_fmt );
155     sprintf( buffer, "%14d%14d%14d%14d%14d%10c", (ptr_lines+ind_lines+val_lines), ptr_lines, ind_lines, val_lines, 0, ' ' );
156     buffer[80] = '\0';
157     fprintf( fp, "%s\n", buffer );
158    
159     /* line 3 */
160     sprintf( buffer, "%c%c%c%11c%14d%14d%14d%14d%10c", 'R', 'U', 'A', ' ', M, N, nz, 0, ' ' );
161     buffer[80] = '\0';
162     fprintf( fp, "%s\n", buffer );
163    
164     /* line 4 */
165     sprintf( buffer, "%16s%16s%20s%28c", ptr_pfmt, ind_pfmt, val_pfmt, ' ');
166     buffer[80]='\0';
167     fprintf( fp, "%s\n", buffer );
168    
169     /* line 5 */
170     /* NOT PRESENT */
171    
172     /* write the actual data */
173     print_data( fp, ptr_perline, ptr_width, (N+1), ptr_fmt, col_ptr, 1, 1 );
174     print_data( fp, ind_perline, ind_width, nz, ind_fmt, row_ind, 1, 1 );
175     print_data( fp, val_perline, val_width, nz, val_fmt, val, 0, 0 );
176     }
177    
178     void Paso_SystemMatrix_saveHB( Paso_SystemMatrix *A_p, char *filename_p )
179     {
180 gross 686 int i, curr_col,j ;
181     int iPtr, iCol, ir, ic;
182 gross 415 index_t index_offset=(A_p->type & MATRIX_FORMAT_OFFSET1 ? 1:0);
183 jgs 150 /* open the file */
184     FILE *fileHandle_p = fopen( filename_p, "w" );
185     if( fileHandle_p == NULL )
186     {
187     Paso_setError(IO_ERROR,"File could not be opened for writing.");
188     return;
189     }
190    
191 gross 415
192 jgs 150 /* check the internal storage type
193     * currently ONLY support for CSC
194     */
195 gross 415 if ( A_p->type & MATRIX_FORMAT_CSC) {
196 jgs 150 if( A_p->row_block_size == 1 && A_p->col_block_size == 1 )
197     {
198     M = A_p->num_rows;
199     N = A_p->num_cols;
200     nz = A_p->len;
201     generate_HB( fileHandle_p, A_p->pattern->ptr, A_p->pattern->index, A_p->val );
202     }
203     else
204     {
205    
206     // fprintf( fileHandle_p, "NEED unrolling!\n" );
207     M = A_p->num_rows*A_p->row_block_size;
208     N = A_p->num_cols*A_p->col_block_size;
209     nz = A_p->len;
210    
211     dim_t *row_ind = MEMALLOC( nz, dim_t );
212     dim_t *col_ind = MEMALLOC( nz, dim_t );
213    
214     i = 0;
215     for( iCol=0; iCol<A_p->pattern->n_ptr; iCol++ )
216     for( ic=0; ic<A_p->col_block_size; ic++)
217 gross 415 for( iPtr=A_p->pattern->ptr[iCol]-index_offset; iPtr<A_p->pattern->ptr[iCol+1]-index_offset; iPtr++)
218 jgs 150 for( ir=0; ir<A_p->row_block_size; ir++ )
219     {
220 gross 415 row_ind[i] = (A_p->pattern->index[iPtr]-index_offset)*A_p->row_block_size+ir+1;
221 jgs 150 col_ind[i] = iCol*A_p->col_block_size+ic+1;
222     i++;
223     }
224    
225     /* get the col_ptr */
226     dim_t *col_ptr = MEMALLOC( (N+1), dim_t );
227    
228     curr_col = 0;
229 gross 686 for(j=0; (j<nz && curr_col<N); curr_col++ )
230 jgs 150 {
231     while( col_ind[j] != curr_col )
232     j++;
233     col_ptr[curr_col] = j;
234     }
235     col_ptr[N] = nz;
236    
237     /* generate the HB file */
238     generate_HB( fileHandle_p, col_ptr, row_ind, A_p->val );
239    
240     /* free the allocated memory */
241     MEMFREE( col_ptr );
242     MEMFREE( col_ind );
243     MEMFREE( row_ind );
244     }
245 gross 415 } else {
246     Paso_setError(TYPE_ERROR,"Paso_SystemMatrix_saveHB: only CSC is currently supported.\n");
247 jgs 150 }
248    
249     /* close the file */
250     fclose( fileHandle_p );
251    
252     return;
253     }
254    
255     /*
256     * $Log$
257     * Revision 1.2 2005/09/15 03:44:39 jgs
258     * Merge of development branch dev-02 back to main trunk on 2005-09-15
259     *
260     * Revision 1.1.2.1 2005/09/05 06:29:48 gross
261     * These files have been extracted from finley to define a stand alone libray for iterative
262     * linear solvers on the ALTIX. main entry through Paso_solve. this version compiles but
263     * has not been tested yet.
264     *
265     *
266     */

Properties

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

  ViewVC Help
Powered by ViewVC 1.1.26