/[escript]/branches/doubleplusgood/paso/src/SparseMatrix_saveHB.cpp
ViewVC logotype

Contents of /branches/doubleplusgood/paso/src/SparseMatrix_saveHB.cpp

Parent Directory Parent Directory | Revision Log Revision Log


Revision 4324 - (show annotations)
Wed Mar 20 00:55:44 2013 UTC (6 years, 1 month ago) by jfenwick
File size: 6404 byte(s)
continues
1
2 /*****************************************************************************
3 *
4 * Copyright (c) 2003-2013 by University of Queensland
5 * http://www.uq.edu.au
6 *
7 * Primary Business: Queensland, Australia
8 * Licensed under the Open Software License version 3.0
9 * http://www.opensource.org/licenses/osl-3.0.php
10 *
11 * Development until 2012 by Earth Systems Science Computational Center (ESSCC)
12 * Development since 2012 by School of Earth Sciences
13 *
14 *****************************************************************************/
15
16
17 /************************************************************************************/
18
19 /* Paso: SparseMatrix saving to Harwell-Boeing format */
20
21 /************************************************************************************/
22
23 /* Copyright: ACcESS Australia 2005 */
24 /* Author: imran@esscc.uq.edu.au */
25
26 /************************************************************************************/
27
28 #include "Paso.h"
29 #include "SparseMatrix.h"
30
31 /* TODO: Refactor the stuff in here into hbio, just like mmio! */
32
33 static dim_t M, N, nz;
34
35 static int calc_digits( int );
36 static void fmt_str( int, int, int*, int*, int*, char*, char* );
37 static void print_data( FILE*, int, int, int, char*, void*, int, int );
38 static void generate_HB( FILE*, dim_t*, dim_t*, double* );
39
40 /* function to get number of digits in an integer */
41 int calc_digits( int var )
42 {
43 int digits = 1;
44 while( (var/=10) )
45 digits++;
46
47 return digits;
48 }
49
50 /* function to generate the format string.
51 *
52 * use maxlen to determine no. of entries per line
53 * use nvalues to determine no. of lines
54 */
55 void fmt_str( int nvalues, int integer, int *width, int *nlines, int *nperline, char *pfmt, char *fmt )
56 {
57 int per_line;
58 int maxlen = *width;
59
60 if( integer && maxlen < 10 )
61 maxlen = 10;
62 else
63 maxlen = 13;
64
65 per_line = 80 / maxlen;
66 *nlines = nvalues / per_line;
67 if( nvalues % per_line )
68 (*nlines)++;
69
70 *nperline = per_line;
71 if( integer )
72 {
73 sprintf( pfmt, "(%dI%d)", per_line, maxlen );
74 sprintf( fmt, "%%%dd", maxlen );
75 }
76 else
77 {
78 sprintf( pfmt, "(1P%dE%d.6)", per_line, maxlen );
79 sprintf( fmt, "%%%d.6E", maxlen );
80 }
81 *width = maxlen;
82 }
83
84 /* function to print the actual data in the right format */
85 void print_data( FILE *fp, int n_perline, int width, int nval, char *fmt, void *ptr, int integer, int adjust )
86 {
87 double *data = reinterpret_cast<double*>(ptr);
88 int entries_done = 0;
89 int padding, i;
90 char pad_fmt[10];
91
92 padding = 80 - n_perline*width;
93 sprintf( pad_fmt, "%%%dc", padding );
94
95 if( adjust != 1 )
96 adjust = 0;
97
98 if( integer )
99 {
100 dim_t *data = reinterpret_cast<dim_t*>(ptr);
101 for(i=0; i<nval; i++ )
102 {
103 fprintf( fp, fmt, data[i]+adjust );
104 entries_done++;
105 if( entries_done == n_perline )
106 {
107 if( padding )
108 fprintf( fp, pad_fmt, ' ' );
109 fprintf( fp, "\n" );
110 entries_done = 0;
111 }
112 }
113 }
114 else
115 {
116 for(i=0; i<nval; i++ )
117 {
118 fprintf( fp, fmt, data[i] );
119 entries_done++;
120 if( entries_done == n_perline )
121 {
122 if( padding )
123 fprintf( fp, pad_fmt, ' ' );
124 fprintf( fp, "\n" );
125 entries_done = 0;
126 }
127 }
128 }
129 if( entries_done )
130 {
131 sprintf( pad_fmt, "%%%dc\n", (80 - entries_done*width) );
132 fprintf( fp, pad_fmt, ' ' );
133 }
134 }
135
136 void generate_HB( FILE *fp, dim_t *col_ptr, dim_t *row_ind, double *val )
137 {
138 char buffer[81];
139
140 int val_lines, ind_lines, ptr_lines;
141 int val_perline, ind_perline, ptr_perline;
142 int val_width, ind_width, ptr_width;
143 char ptr_pfmt[7], ind_pfmt[7], val_pfmt[11];
144 char ptr_fmt[10], ind_fmt[10], val_fmt[10];
145
146 /* line 1 */
147 sprintf( buffer, "%-72s%-8s", "Matrix Title", "Key" );
148 buffer[80] = '\0';
149 fprintf( fp, "%s\n", buffer );
150
151 /* line 2 */
152 ptr_width = calc_digits( nz+1 );
153 fmt_str( N+1, 1, &ptr_width, &ptr_lines, &ptr_perline, ptr_pfmt, ptr_fmt );
154 ind_width = calc_digits( N );
155 fmt_str( nz, 1, &ind_width, &ind_lines, &ind_perline, ind_pfmt, ind_fmt );
156 val_width = 13;
157 fmt_str( nz, 0, &val_width, &val_lines, &val_perline, val_pfmt, val_fmt );
158 sprintf( buffer, "%14d%14d%14d%14d%14d%10c", (ptr_lines+ind_lines+val_lines), ptr_lines, ind_lines, val_lines, 0, ' ' );
159 buffer[80] = '\0';
160 fprintf( fp, "%s\n", buffer );
161
162 /* line 3 */
163 sprintf( buffer, "%c%c%c%11c%14d%14d%14d%14d%10c", 'R', 'U', 'A', ' ', M, N, nz, 0, ' ' );
164 buffer[80] = '\0';
165 fprintf( fp, "%s\n", buffer );
166
167 /* line 4 */
168 sprintf( buffer, "%16s%16s%20s%28c", ptr_pfmt, ind_pfmt, val_pfmt, ' ');
169 buffer[80]='\0';
170 fprintf( fp, "%s\n", buffer );
171
172 /* line 5 */
173 /* NOT PRESENT */
174
175 /* write the actual data */
176 print_data( fp, ptr_perline, ptr_width, (N+1), ptr_fmt, col_ptr, 1, 1 );
177 print_data( fp, ind_perline, ind_width, nz, ind_fmt, row_ind, 1, 1 );
178 print_data( fp, val_perline, val_width, nz, val_fmt, val, 0, 0 );
179 }
180
181 void Paso_SparseMatrix_saveHB_CSC( Paso_SparseMatrix *A_p, FILE* fileHandle_p ) {
182 int i, curr_col,j ;
183 int iPtr, iCol, ir, ic;
184 index_t index_offset=(A_p->type & MATRIX_FORMAT_OFFSET1 ? 1:0);
185 dim_t *row_ind=NULL, *col_ind = NULL, *col_ptr=NULL;
186 int nz = A_p->len;
187
188 if (A_p->val == NULL) {
189 Esys_setError(TYPE_ERROR,"Paso_SparseMatrix_saveHB_CSC: unsupported format detected.\n");
190 return;
191 }
192
193 if( A_p->row_block_size == 1 && A_p->col_block_size == 1 ) {
194 M = A_p->numRows;
195 N = A_p->numCols;
196 generate_HB( fileHandle_p, A_p->pattern->ptr, A_p->pattern->index, A_p->val );
197 } else {
198 M = A_p->numRows*A_p->row_block_size;
199 N = A_p->numCols*A_p->col_block_size;
200
201 row_ind = new dim_t [nz];
202 col_ind = new dim_t [nz];
203
204 i = 0;
205 for( iCol=0; iCol<A_p->pattern->numOutput; iCol++ )
206 for( ic=0; ic<A_p->col_block_size; ic++)
207 for( iPtr=A_p->pattern->ptr[iCol]-index_offset; iPtr<A_p->pattern->ptr[iCol+1]-index_offset; iPtr++)
208 for( ir=0; ir<A_p->row_block_size; ir++ ) {
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 /* get the col_ptr */
214 col_ptr = new dim_t [(N+1)];
215
216 curr_col = 0;
217 for(j=0; (j<nz && curr_col<N); curr_col++ ) {
218 while( col_ind[j] != curr_col )
219 j++;
220 col_ptr[curr_col] = j;
221 }
222 col_ptr[N] = nz;
223
224 /* generate the HB file */
225 generate_HB( fileHandle_p, col_ptr, row_ind, A_p->val );
226
227 /* free the allocated memory */
228 delete[] col_ptr;
229 delete[] col_ind;
230 delete[] row_ind;
231 }
232 return;
233 }

  ViewVC Help
Powered by ViewVC 1.1.26