/[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 150 - (hide annotations)
Thu Sep 15 03:44:45 2005 UTC (14 years, 1 month ago) by jgs
Original Path: trunk/esys2/paso/src/SystemMatrix_saveHB.c
File MIME type: text/plain
File size: 6597 byte(s)
Merge of development branch dev-02 back to main trunk on 2005-09-15

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

Properties

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

  ViewVC Help
Powered by ViewVC 1.1.26