Hitmap 1.3
 All Data Structures Namespaces Files Functions Variables Typedefs Friends Macros Groups Pages
hit_file.c
Go to the documentation of this file.
1 
11 /*
12  * <license>
13  *
14  * Hitmap v1.2
15  *
16  * This software is provided to enhance knowledge and encourage progress in the scientific
17  * community. It should be used only for research and educational purposes. Any reproduction
18  * or use for commercial purpose, public redistribution, in source or binary forms, with or
19  * without modifications, is NOT ALLOWED without the previous authorization of the copyright
20  * holder. The origin of this software must not be misrepresented; you must not claim that you
21  * wrote the original software. If you use this software for any purpose (e.g. publication),
22  * a reference to the software package and the authors must be included.
23  *
24  * THIS SOFTWARE IS PROVIDED BY THE COPYRIGHT HOLDER AND CONTRIBUTORS "AS IS" AND ANY
25  * EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED TO, THE IMPLIED WARRANTIES OF
26  * MERCHANTABILITY AND FITNESS FOR A PARTICULAR PURPOSE ARE DISCLAIMED. IN NO EVENT SHALL
27  * THE AUTHORS OR CONTRIBUTORS BE LIABLE FOR ANY DIRECT, INDIRECT, INCIDENTAL, SPECIAL,
28  * EXEMPLARY, OR CONSEQUENTIAL DAMAGES (INCLUDING, BUT NOT LIMITED TO, PROCUREMENT OF
29  * SUBSTITUTE GOODS OR SERVICES; LOSS OF USE, DATA, OR PROFITS; OR BUSINESS INTERRUPTION)
30  * HOWEVER CAUSED AND ON ANY THEORY OF LIABILITY, WHETHER IN CONTRACT, STRICT LIABILITY,
31  * OR TORT (INCLUDING NEGLIGENCE OR OTHERWISE) ARISING IN ANY WAY OUT OF THE USE OF THIS
32  * SOFTWARE, EVEN IF ADVISED OF THE POSSIBILITY OF SUCH DAMAGE.
33  *
34  * Copyright (c) 2007-2015, Trasgo Group, Universidad de Valladolid.
35  * All rights reserved.
36  *
37  * More information on http://trasgo.infor.uva.es/
38  *
39  * </license>
40 */
41 
42 #include <hit_file.h>
43 #include <stdlib.h>
44 #include <stdio.h>
45 #include <errno.h>
46 #include <ctype.h>
47 
48 #include <hit_error.h>
49 #include <hit_allocP.h>
50 #include <hit_funcop.h>
51 #include <hit_cshape.h>
52 #include <hit_bshape.h>
53 
54 #ifndef INCLUDE_MATRIX_IO__H
55  #define INCLUDE_MATRIX_IO__H
56  #include "matrix_io.h"
57 #endif
58 
66 int compare_idxtype (const void * a, const void * b){
67  return (*(const idxtype *)a - *(const idxtype *)b);
68 }
69 
70 
72 void hit_csr_renameIndexes(int Nr, int * xadj, int * adjncy){
73 
74  int i;
75  for(i=0;i<=Nr;i++){
76  xadj[i]--;
77  }
78 
79  for(i=0;i<xadj[Nr];i++){
80  adjncy[i]--;
81  }
82 }
83 
84 
86 int hit_csr_renameIndexes_and_Symmetrize(int * pNr, int * pNc, int nz, int ** pxadj, int ** padjncy,
87  const char * file, int line){
88 
89  // 0. Original CSR
90  int i,j;
91  int *colptr = *pxadj;
92  int *rowind = *padjncy;
93 
94  // 1. Get the number of vertices
95  int N = hit_max(*pNr,*pNc);
96 
97  //printf("[%d,%d]: max: %d\n",*pNr,*pNc,N);
98 
99  if(*pNr < N){
100 
101  // @arturo Ago 2015: New allocP interface
102  // hit_realloc(colptr,(size_t) (N+1) * sizeof(idxtype),idxtype*);
103  hit_realloc(colptr, idxtype, N+1 );
104 
105  if (colptr == NULL) {
106  hit_errInternal(__FUNCTION__,"Symmetrize: not enough memory for data structures","",file,line);
107  return HIT_ERROR;
108  }
109  for(i=(*pNr)+1;i<N+1;i++){
110  colptr[i] = colptr[i-1];
111  }
112  }
113  *pNr = N;
114  *pNc = N;
115 
116 
117 
118 
119  // 2. Allocate memory for CSR format
120  (*pxadj) = idxmalloc(N+1, (char *)&"Symmetrize: xadj");
121  (*padjncy) = idxmalloc(2*nz, (char *)&"Symmetrize: adjncy");
122 
123  // 3. Allocate memory for edge counters
124  int * n_edge;
125  // @arturo Ago 2015: New allocP interface
126  // hit_calloc(n_edge,(size_t)N,sizeof(int),int*);
127  hit_calloc(n_edge, int, N );
128 
129  // 4. Check allocations.
130  if ((*pxadj)==NULL || (*padjncy)==NULL || n_edge==NULL) {
131  hit_errInternal(__FUNCTION__,"Symmetrize: not enough memory for data structures","",file,line);
132  return HIT_ERROR;
133  }
134 
135  // 5. Count the number of edges that do not appear in the HB (unsymmetric)
136  for(i=0;i<N;i++){
137 
138  int begin = colptr[i]-1;
139  int end = colptr[i+1]-1;
140  for(j=begin;j<end;j++){
141 
142  int target = rowind[j]-1;
143  int begin2 = colptr[target]-1;
144  int end2 = colptr[target+1]-1;
145  int j2;
146  int notfound = 1;
147 
148  for(j2=begin2;j2<end2;j2++){
149  int target2 = rowind[j2]-1;
150  if(target2 == i){
151  notfound = 0;
152  break;
153  }
154  }
155  n_edge[target] += notfound;
156  }
157  }
158 
159 
160  // 6. Init the pxadj with the correct values
161  (*pxadj)[0] = 0;
162  for(i=0;i<N;i++){
163  int begin = colptr[i]-1;
164  int end = colptr[i+1]-1;
165  int length = (end-begin) + n_edge[i];
166 
167  (*pxadj)[i+1] = (*pxadj)[i] + length;
168  }
169 
170  // 7. Copy the HB unsymmetric edges in adjncy
171  for(i=0;i<N;i++){
172 
173  int begin = colptr[i]-1;
174  int end = colptr[i+1]-1;
175  int length = (end-begin);
176 
177  int newbegin = (*pxadj)[i];
178 
179  for(j=0;j<length;j++){
180 
181  int target = rowind[begin+j]-1;
182  (*padjncy)[newbegin+j] = target;
183  }
184  }
185 
186  // 8. Copy the remain edges
187  for(i=0;i<N;i++){
188 
189  int begin = colptr[i]-1;
190  int end = colptr[i+1]-1;
191  for(j=begin;j<end;j++){
192 
193  int target = rowind[j]-1;
194  int begin2 = colptr[target]-1;
195  int end2 = colptr[target+1]-1;
196  int j2;
197  int notfound = 1;
198 
199  // Optimization to skip complete edges
200  if(n_edge[target] == 0) continue;
201 
202  for(j2=begin2;j2<end2;j2++){
203  int target2 = rowind[j2]-1;
204  if(target2 == i){
205  notfound = 0;
206  break;
207  }
208  }
209 
210  // Copy the edge
211  if(notfound){
212  int pos = (*pxadj)[target+1]-n_edge[target];
213  (*padjncy)[pos] = i;
214  n_edge[target]--;
215  }
216  }
217  }
218 
219 
220  // 9. Reorder the edges
221  for(i=0;i<N;i++){
222 
223  int begin = (*pxadj)[i];
224  int end = (*pxadj)[i+1];
225 
226  qsort(&(*padjncy)[begin], (size_t)(end-begin), sizeof(idxtype), compare_idxtype);
227  }
228 
229  // 10. Free extra memory
230  // @arturo Ago 2015: New allocP interface
231  // hit_realloc(*padjncy,(size_t) ((*pxadj)[N]) * sizeof(idxtype),idxtype*);
232  hit_realloc(*padjncy, idxtype, (*pxadj)[N] );
233 
234  // 99. FREE DATA STRUCTURES
235  free(colptr);
236  free(rowind);
237  free(n_edge);
238 
239 
240  return HIT_OK;
241 
242 }
243 
244 
245 
246 HitShape hit_fileMMRead_toCSR_Internal(const char *fileName, int create_graph, int rank,
247  const char * file, int line){
248 
249  // 0. If the current processor is not root return a null shape.
250  if(rank != 0){
251  hit_warnInternal(__FUNCTION__,"Only root processor should use this function","",file,line);
252  return HIT_CSR_SHAPE_NULL;
253  }
254 
255  // 1. READ HB HEADER
256  // Rows and columns sizes
257  int Nr, Nc;
258  // Nonzero elements
259  int nz;
260  // Matrix type
261  char * type;
262  // Number of right-hand-sides
263  int Nrhs;
264 
265  // 2. READ THE MATRIX
266 
267  // Row and column pointer.
268  int *xadj = NULL, *adjncy = NULL;
269  // Array with the values.
270  double *val = NULL;
271 
272  // 2.1 Read the matrix allocating storage arrays
273 
274  read_mm_matrix ((char *)fileName, &Nr, &Nc, &nz, &xadj, &adjncy, &val);
275 
276 
277  // 3. TRANSFORM THE CSR TO 0-index AND SIMETRIZE
278  if( create_graph){
279  // hit_csr_renameIndexes_and_Symmetrize(&Nr, &Nc, nz, &xadj, &adjncy, file, line);
280  } else {
281  //hit_csr_renameIndexes(Nr, xadj, adjncy);
282  }
283 
284 
285 #ifdef DEBUG
286  {
287  int i,j;
288  printf("\nMM: Nr: %d Nc: %d \n", Nr, Nc);
289  printf("\t");
290  for(i=0;i<Nr;i++){
291  int begin = xadj[i];
292  int end = xadj[i+1];
293  for(j=begin;j<end;j++){
294  printf(" %d ", adjncy[j]);
295  }
296  printf("|\n\t");
297  }
298  printf("\n");
299  }
300 #endif
301 
302 
303  // 4. CREATE THE NAME LISTS
304  // Row and Column lists.
305  HitNameList listR;
306  HitNameList listC;
307 
308  hit_nameListCreate(&listR, Nr);
309 
310  if(create_graph){
311  listC = listR;
312  } else {
313  hit_nameListCreate(&listC, Nc);
314  }
315 
316 
317  // 5. BUILD THE SHAPE
319  hit_cShapeCard(shape,0) = Nr;
320  hit_cShapeCard(shape,1) = Nc;
321  hit_cShapeXadj(shape) = xadj;
322  hit_cShapeAdjncy(shape) = adjncy;
323  hit_cShapeNameList(shape,0) = listR;
324  hit_cShapeNameList(shape,1) = listC;
325 
326 
327  // 99. Free resources
328  free(type);
329  free(val);
330 
331  return shape;
332 }
333 
334 
335 
336 
337 
338 
339 
340 HitShape hit_fileHBRead_toCSR_Internal(const char *fileName, int create_graph, int rank,
341  const char * file, int line){
342 
343  // 0. If the current processor is not root return a null shape.
344  if(rank != 0){
345  hit_warnInternal(__FUNCTION__,"Only root processor should use this function","",file,line);
346  return HIT_CSR_SHAPE_NULL;
347  }
348 
349  // 1. READ HB HEADER
350  // Rows and columns sizes
351  int Nr, Nc;
352  // Nonzero elements
353  int nz;
354  // Matrix type
355  char * type;
356  // Number of right-hand-sides
357  int Nrhs;
358 
359  // 1.1 Read HB header
360  if( !readHB_info(fileName, &Nr, &Nc, &nz, &type, &Nrhs) ) {
361  hit_errInternal(__FUNCTION__,"Reading Harwell-Boeing file header","",file,line);
362  return HIT_CSR_SHAPE_NULL;
363  }
364 
365 #ifdef DEBUG
366  fprintf(stderr,"[HB file] Nr: %d, Nc: %d, nz: %d, type: %s, Nrhs: %d\n",Nr,Nc,nz,type,Nrhs);
367 #endif
368 
369 
370  // 1.2. Reject Righ-hand side matrices
371  if( Nrhs > 0 ) {
372  hit_warnInternal(__FUNCTION__,"Hitmap ignores Harwell-Boeing matrices with right-hand-sides","",file,line);
373  }
374 
375  // 1.3 Reject not pattern matrices
376  // (P)attern, (R)eal or (C)omplex.
377  if( type[0] != 'P' ){
378  hit_warnInternal(__FUNCTION__,"Hitmap ignores Harwell-Boeing matrix values","",file,line);
379  }
380 
381  // 1.4 Reject not assembled matrices
382  if( type[1] != 'U' && type[1] != 'R' && type[1] != 'S' ){
383  hit_errInternal(__FUNCTION__,"Harwell-Boeing matrix type not supported","",file,line);
384  return HIT_CSR_SHAPE_NULL;
385  }
386 
387  // 1.5 Reject not assembled matrices
388  if( type[2] != 'A' ){
389  hit_errInternal(__FUNCTION__,"Harwell-Boeing matrix is not assembled","",file,line);
390  return HIT_CSR_SHAPE_NULL;
391  }
392 
393 
394 
395  // 2. READ THE MATRIX
396 
397  // Row and column pointer.
398  int *xadj = NULL, *adjncy = NULL;
399  // Array with the values.
400  double *val = NULL;
401 
402  // 2.1 Read the matrix allocating storage arrays
403  if( !readHB_newmat_double(fileName, &Nr, &Nc, &nz, &xadj, &adjncy, &val) ) {
404  hit_errInternal(__FUNCTION__,"Reading Harwell-Boeing matrix data","",file,line);
405  return HIT_CSR_SHAPE_NULL;
406  }
407 
408 
409  // 3. TRANSFORM THE CSR TO 0-index AND SIMETRIZE
410  if(type[1] == 'S' || create_graph){
411  hit_csr_renameIndexes_and_Symmetrize(&Nr, &Nc, nz, &xadj, &adjncy, file, line);
412  } else {
413  hit_csr_renameIndexes(Nr, xadj, adjncy);
414  }
415 
416 #ifdef DEBUG
417  {
418  int i,j;
419  printf("\nHB: Nr: %d Nc: %d \n", Nr, Nc);
420  printf("\t");
421  for(i=0;i<Nr;i++){
422  int begin = xadj[i];
423  int end = xadj[i+1];
424  for(j=begin;j<end;j++){
425  printf(" %d ", adjncy[j]);
426  }
427  printf("|\n\t");
428  }
429  printf("\n");
430  }
431 #endif
432 
433  // 4. CREATE THE NAME LISTS
434  // Row and Column lists.
435  HitNameList listR;
436  HitNameList listC;
437 
438  hit_nameListCreate(&listR, Nr);
439 
440  if(create_graph){
441  listC = listR;
442  } else {
443  hit_nameListCreate(&listC, Nc);
444  }
445 
446 
447  // 5. BUILD THE SHAPE
449  hit_cShapeCard(shape,0) = Nr;
450  hit_cShapeCard(shape,1) = Nc;
451  hit_cShapeXadj(shape) = xadj;
452  hit_cShapeAdjncy(shape) = adjncy;
453  hit_cShapeNameList(shape,0) = listR;
454  hit_cShapeNameList(shape,1) = listC;
455 
456 
457  // 99. Free resources
458  free(type);
459  free(val);
460 
461  return shape;
462 }
463 
464 
465 
466 HitShape hit_fileHBRead_toBitmap_Internal(const char *fileName, int create_graph, int rank,
467  const char * file, int line){
468 
469  // 0. If the current processor is not root return a null shape.
470  if(rank != 0){
471  hit_warnInternal(__FUNCTION__,"Only root processor should use this function","",file,line);
472  return HIT_CSR_SHAPE_NULL;
473  }
474 
475  // 1. READ HB HEADER
476  // Rows and columns sizes
477  int Nr, Nc;
478  // Nonzero elements
479  int nz;
480  // Matrix type
481  char * type;
482  // Number of right-hand-sides
483  int Nrhs;
484 
485  // 1.1 Read HB header
486  if( !readHB_info(fileName, &Nr, &Nc, &nz, &type, &Nrhs) ) {
487  hit_errInternal(__FUNCTION__,"Reading Harwell-Boeing file header","",file,line);
488  return HIT_CSR_SHAPE_NULL;
489  }
490 
491 #ifdef DEBUG
492  fprintf(stderr,"[HB file] Nr: %d, Nc: %d, nz: %d, type: %s, Nrhs: %d\n",Nr,Nc,nz,type,Nrhs);
493 #endif
494 
495 
496  // 1.2. Reject Righ-hand side matrices
497  if( Nrhs > 0 ) {
498  hit_warnInternal(__FUNCTION__,"Hitmap ignores Harwell-Boeing matrices with right-hand-sides","",file,line);
499  }
500 
501  // 1.3 Reject not pattern matrices
502  // (P)attern, (R)eal or (C)omplex.
503  if( type[0] != 'P' ){
504  hit_warnInternal(__FUNCTION__,"Hitmap ignores Harwell-Boeing matrix values","",file,line);
505  }
506 
507  // 1.4 Reject not assembled matrices
508  if( type[1] != 'U' && type[1] != 'R' && type[1] != 'S' ){
509  hit_errInternal(__FUNCTION__,"Harwell-Boeing matrix type not supported","",file,line);
510  return HIT_CSR_SHAPE_NULL;
511  }
512 
513  // 1.5 Reject not assembled matrices
514  if( type[2] != 'A' ){
515  hit_errInternal(__FUNCTION__,"Harwell-Boeing matrix is not assembled","",file,line);
516  return HIT_CSR_SHAPE_NULL;
517  }
518 
519 
520 
521  // 2. READ THE MATRIX
522 
523  // Row and column pointer.
524  int *xadj = NULL, *adjncy = NULL;
525  // Array with the values.
526  double *val = NULL;
527 
528  // 2.1 Read the matrix allocating storage arrays
529  if( !readHB_newmat_double(fileName, &Nr, &Nc, &nz, &xadj, &adjncy, &val) ) {
530  hit_errInternal(__FUNCTION__,"Reading Harwell-Boeing matrix data","",file,line);
531  return HIT_CSR_SHAPE_NULL;
532  }
533 
534 
535  /* 3. LOAD THE MATRIX IN THE SHAPE */
536 
537  /* 3.1 New HitBShape */
539 
540  if(create_graph){
541  shape = hit_bitmapShape(Nr);
542  } else {
543  shape = hit_bitmapShapeMatrix(Nr,Nc);
544  }
545 
546  /* 3.2 Fill the tile */
547  int i,j;
548  for(i=0; i<Nr; i++) {
549  int cstyleIni = xadj[i]-1;
550  int cstyleEnd = xadj[i+1]-1;
551 
552  for(j=cstyleIni; j<cstyleEnd; j++) {
553  hit_bShapeSet(shape,i,adjncy[j]-1);
554  if(type[1] == 'S' || create_graph){
555  hit_bShapeSet(shape,adjncy[j]-1,i);
556  }
557  }
558  }
559 
560  // 99. Free resources
561  free(type);
562  free(val);
563  free(xadj);
564  free(adjncy);
565 
566  return shape;
567 }
568 
569 
570 
571 
572 
573 void hit_fileHBWriteInternal(const char * hbfile, HitShape shape, int rank, const char * cfile, int cline) {
574  HIT_NOT_USED( rank );
575 
576 // 10 x 8 = 80 characters per line.
577 #define NUMBERS_LINE 10
578 #define NUMBERS_SIZE 8
579 
580 
581  if(hit_shapeType(shape) != HIT_CSR_SHAPE){
582  hit_errInternal(__FUNCTION__,"Wrong shape type","",cfile,cline);
583  }
584 
585 
586  FILE * file = fopen(hbfile, "w");
587 
588  // Number of Vertices and edges
589  int nvertices = hit_cShapeNvertices(shape);
590  int nedges = hit_cShapeNedges(shape);
591 
592  // LINE 1
593  // 1.1 Title (72 characters)
594  fprintf(file,"%-72s","HitMap HB Export");
595  // 1.2, (8 characters)
596  fprintf(file,"%-8s","key");
597  fprintf(file,"\n");
598 
599  // LINE 2
600  //2.1 TOTCRD, total number of data lines, (14 characters)
601  //2.2 PTRCRD, number of data lines for pointers, (14 characters)
602  //2.3 INDCRD, number of data lines for row or variable indices, (14 characters)
603  int ptrcrd = (int) ceil(((double) nvertices + 1) / NUMBERS_LINE );
604  int indcrd = (int) ceil((double) nedges / NUMBERS_LINE );
605  fprintf(file,"%14d",ptrcrd + indcrd);
606  fprintf(file,"%14d",ptrcrd);
607  fprintf(file,"%14d",indcrd);
608  //2.4 VALCRD, number of data lines for numerical values of matrix entries, (14 characters)
609  //2.5 RHSCRD, number of data lines for right hand side vectors, starting guesses, and solutions, (14 characters)
610  fprintf(file,"%14d",0);
611  fprintf(file,"%14s","");
612  fprintf(file,"\n");
613 
614  // LINE 2
615  // MXTYPE, matrix type, (3 characters)
616  // pua: Pattern Unsymmetric Assembled
617  fprintf(file,"%3s","pua");
618  // blank space, (11 characters)
619  fprintf(file,"%11s","");
620  // NROW, integer, number of rows or variables, (14 characters)
621  fprintf(file,"%14d", nvertices);
622  // NCOL, integer, number of columns or elements, (14 characters)
623  fprintf(file,"%14d", nvertices);
624  // NNZERO, number of nonzero entries. (14 characters)
625  fprintf(file,"%14d", nedges);
626  // NELTVL, this is 0. (14 characters)
627  fprintf(file,"%14d", 0);
628  fprintf(file,"\n");
629 
630 
631  // LINE 4
632  // PTRFMT, FORTRAN I/O format for pointers, (16 characters)
633  fprintf(file,"%-16s","(10I8)");
634  // INDFMT, FORTRAN I/O format for row or variable indices, (16 characters)
635  fprintf(file,"%-16s","(10I8)");
636  // VALFMT, FORTRAN I/O format for matrix entries, (20 characters)
637  fprintf(file,"%-20s","");
638  // RHSFMT, FORTRAN I/O format for right hand sides, initial guesses, and solutions, (20 characters)
639  fprintf(file,"%-20s","");
640  //fprintf(file,"\n");
641 
642  int i;
643 
644  for(i=0;i<=nvertices;i++){
645 
646  if(i % NUMBERS_LINE == 0 ) fprintf(file,"\n");
647  int vertex = hit_cShapeXadj(shape)[i];
648  fprintf(file, "%8d",vertex+1);
649  }
650 
651 
652  for(i=0;i<nedges;i++){
653 
654  if(i % NUMBERS_LINE == 0 ) fprintf(file,"\n");
655  int edge = hit_cShapeAdjncy(shape)[i];
656  fprintf(file, "%8d",edge+1);
657  }
658 
659  fclose(file);
660 
661 }
662 
663 
664 void hit_fileHBWriteBitmapInternal(const char * hbfile, HitShape shape, int rank, const char * cfile, int cline){
665  HIT_NOT_USED( rank );
666 
667 // 10 x 8 = 80 characters per line.
668 #define NUMBERS_LINE 10
669 #define NUMBERS_SIZE 8
670 
671 
672  if(hit_shapeType(shape) != HIT_BITMAP_SHAPE){
673  hit_errInternal(__FUNCTION__,"Wrong shape type","",cfile,cline);
674  }
675 
676 
677  FILE * file = fopen(hbfile, "w");
678  //FILE * file = stdout;
679 
680  // Number of Vertices and edges
681  int nvertices = hit_bShapeNvertices(shape);
682  int nedges = hit_bShapeNedges(shape);
683 
684  // LINE 1
685  // 1.1 Title (72 characters)
686  fprintf(file,"%-72s","HitMap HB Export");
687  // 1.2, (8 characters)
688  fprintf(file,"%-8s","key");
689  fprintf(file,"\n");
690 
691  // LINE 2
692  //2.1 TOTCRD, total number of data lines, (14 characters)
693  //2.2 PTRCRD, number of data lines for pointers, (14 characters)
694  //2.3 INDCRD, number of data lines for row or variable indices, (14 characters)
695  int ptrcrd = (int) ceil(((double) nvertices + 1) / NUMBERS_LINE );
696  int indcrd = (int) ceil((double) nedges / NUMBERS_LINE );
697  fprintf(file,"%14d",ptrcrd + indcrd);
698  fprintf(file,"%14d",ptrcrd);
699  fprintf(file,"%14d",indcrd);
700  //2.4 VALCRD, number of data lines for numerical values of matrix entries, (14 characters)
701  //2.5 RHSCRD, number of data lines for right hand side vectors, starting guesses, and solutions, (14 characters)
702  fprintf(file,"%14d",0);
703  fprintf(file,"%14s","");
704  fprintf(file,"\n");
705 
706  // LINE 2
707  // MXTYPE, matrix type, (3 characters)
708  // pua: Pattern Unsymmetric Assembled
709  fprintf(file,"%3s","pua");
710  // blank space, (11 characters)
711  fprintf(file,"%11s","");
712  // NROW, integer, number of rows or variables, (14 characters)
713  fprintf(file,"%14d", nvertices);
714  // NCOL, integer, number of columns or elements, (14 characters)
715  fprintf(file,"%14d", nvertices);
716  // NNZERO, number of nonzero entries. (14 characters)
717  fprintf(file,"%14d", nedges);
718  // NELTVL, this is 0. (14 characters)
719  fprintf(file,"%14d", 0);
720  fprintf(file,"\n");
721 
722 
723  // LINE 4
724  // PTRFMT, FORTRAN I/O format for pointers, (16 characters)
725  fprintf(file,"%-16s","(10I8)");
726  // INDFMT, FORTRAN I/O format for row or variable indices, (16 characters)
727  fprintf(file,"%-16s","(10I8)");
728  // VALFMT, FORTRAN I/O format for matrix entries, (20 characters)
729  fprintf(file,"%-20s","");
730  // RHSFMT, FORTRAN I/O format for right hand sides, initial guesses, and solutions, (20 characters)
731  fprintf(file,"%-20s","");
732  //fprintf(file,"\n");
733 
734 
735  int indexId = 0;
736  int vertex;
737  hit_bShapeVertexIterator(vertex,shape){
738 
739  if(vertex % NUMBERS_LINE == 0 ) fprintf(file,"\n");
740  fprintf(file, "%8d",indexId+1);
741 
742  int edge;
743  hit_bShapeEdgeIterator(edge,shape,vertex){
744  indexId++;
745  }
746 
747  }
748  if(vertex % NUMBERS_LINE == 0 ) fprintf(file,"\n");
749  fprintf(file, "%8d",indexId+1);
750 
751  int edge_count = 0;
752  hit_bShapeVertexIterator(vertex,shape){
753 
754  int edge;
755  hit_bShapeEdgeIterator(edge,shape,vertex){
756 
757  if(edge_count % NUMBERS_LINE == 0 ) fprintf(file,"\n");
758  edge_count++;
759 
760  int target = hit_bShapeEdgeTarget(shape,edge);
761 
762  fprintf(file, "%8d",target+1);
763 
764  }
765 
766  }
767 
768 
769 
770  fclose(file);
771 
772 }
773 
774 
775 
776 
777 
778 int hit_fileHBVerticesInternal(const char * hbfile, int rank, const char * file, int line){
779 
780  // Check if the current processor is not root.
781  if(rank != 0){
782  hit_warnInternal(__FUNCTION__,"Only root processor should use this function","",file,line);
783  return 0;
784  }
785 
786  /* 1. READ HB HEADER */
787  // Rows and columns sizes.
788  int M, N;
789  // Nonzero elements,
790  int nz;
791  // Matrix type.
792  char * type;
793  // Number of right-hand-sides
794  int Nrhs;
795 
796  /* 1.1 Read HB header */
797  if( !readHB_info(hbfile, &M, &N, &nz, &type, &Nrhs) ) {
798  hit_errInternal(__FUNCTION__,"Reading Harwell-Boeing file header","",file,line);
799  return HIT_ERROR;
800  }
801 
802  /* 1.3. Reject non-squeare matrices */
803  if( M != N ) {
804  hit_errInternal(__FUNCTION__,"Recognized Harwell-Boeing matrices must be square","",file,line);
805  return HIT_ERROR;
806  }
807 
808  // Free type string.
809  free(type);
810 
811  return N;
812 }
813 
814 
815 
816 int hit_fileHBReadDenseInternal(const char * fileName, int rank, void * tileP, const char * file, int line){
817 
818  // Check if the current processor is not root.
819  if(rank != 0){
820  hit_warnInternal(__FUNCTION__,"Only root processor should use this function","",file,line);
821  return HIT_ERROR;
822  }
823 
824  int i,j;
825 
826  /* 0. CHECK TILE */
827  hit_tileNewType(int);
828 
829  HitTile_int *tile = (HitTile_int *)tileP;
830  if(hit_tileIsNull(*tile)){
831  hit_errInternal(__FUNCTION__,"The tile must be allocated","",file,line);
832  }
833 
834 
835  /* 1. READ HB HEADER */
836  // Rows and columns sizes.
837  int M, N;
838  // Nonzero elements,
839  int nz;
840  // Matrix type.
841  char * type;
842  // Number of right-hand-sides
843  int Nrhs;
844 
845  /* 1.1 Read HB header */
846  if( !readHB_info(fileName, &M, &N, &nz, &type, &Nrhs) ) {
847  hit_errInternal(__FUNCTION__,"Reading Harwell-Boeing file header","",file,line);
848  return HIT_ERROR;
849  }
850 
851 #ifdef DEBUG
852  fprintf(stderr,"DATOS: M: %d, N: %d, nz: %d, type: %s, Nrhs: %d\n",M,N,nz,type,Nrhs);
853 #endif
854 
855  /* 1.2. Reject Righ-hand side matrices */
856  if( Nrhs > 0 ) {
857  hit_errInternal(__FUNCTION__,"Recognized Harwell-Boeing files must not have RHS","",file,line);
858  return HIT_ERROR;
859  }
860 
861  /* 1.3. Reject non-square matrices */
862  if( M != N ) {
863  hit_errInternal(__FUNCTION__,"Recognized Harwell-Boeing matrices must be square","",file,line);
864  return HIT_ERROR;
865  }
866 
867  // Free type string.
868  free(type);
869 
870 
871 
872  /* 2. READ THE MATRIX */
873 
874  // Row and column pointer.
875  int *colptr = NULL, *rowind = NULL;
876  // Array with the values.
877  double *val = NULL;
878 
879 
880  /* 2. Read the matrix allocating storage arrays */
881  if( !readHB_newmat_double(fileName, &M, &N, &nz, &colptr, &rowind, &val) ) {
882  hit_errInternal(__FUNCTION__,"Reading Harwell-Boeing matrix data","",file,line);
883  return HIT_ERROR;
884  }
885 
886 #ifdef DEBUG_LEVEL2
887  /*
888  {int i;
889  for (i=0; i<N+1; i++) {
890  fprintf(stderr,"colptr[%i] = %i\n",i,colptr[i]-1);
891  }
892  for (i=0; i<nz; i++) {
893  fprintf(stderr,"rowind[%i] = %i\n",i,rowind[i]-1);
894  }}
895  */
896  {int i,j;
897  printf("\t");
898  for(i=0;i<N;i++){
899  int begin = colptr[i]-1;
900  int end = colptr[i+1]-1;
901  for(j=begin;j<end;j++){
902  printf(" %d ", rowind[j]-1);
903  }
904  printf("|\n\t");
905  }
906  printf("\n");}
907 #endif
908 
909 
910  /* 3. LOAD THE MATRIX IN THE TILE STRUCTURE */
911 
912  /* 3.1 Clean the tile */
913  int zero = 0;
914  hit_tileFill(tileP,&zero);
915 
916  /* 3.2 Fill the tile */
917  for(i=0; i<N; i++) {
918  int cstyleIni = colptr[i]-1;
919  int cstyleEnd = colptr[i+1]-1;
920 
921  for(j=cstyleIni; j<cstyleEnd; j++) {
922 
923  hit_tileElemAt2(*tile,i,rowind[j]-1) = 1;
924  hit_tileElemAt2(*tile,rowind[j]-1,i) = 1;
925  //Mat[i][rowind[j]-1] = 1;
926  //Mat[rowind[j]-1][i] = 1;
927  }
928  }
929 
930 
931  /* 4. FREE DATA STRUCTURES */
932  free(colptr);
933  free(rowind);
934  free(val);
935 
936  /* END */
937  return HIT_OK;
938 
939 }
940 
941 
942 HitShape hit_fileCSRReadInternal(const char * csrfile, int rank, const char * cfile, int line){
943 
944  // If the current processor is not root return a null shape.
945  if(rank != 0){
946  hit_warnInternal(__FUNCTION__,"Only root processor should use this function","",cfile,line);
947  return HIT_CSR_SHAPE_NULL;
948  }
949 
950  FILE * file = fopen(csrfile,"r");
951  if(file == NULL) return HIT_CSR_SHAPE_NULL;
952 
953 
955  int n;
956 
957  char c = (char) getc(file);
958  while( c != EOF ){
959 
960  if(c == '#' || c == ' ' || c == '\n'){
961  while( c != '\n' ) c = (char) getc(file);
962  } else {
963  ungetc(c,file);
964  break;
965  }
966  c = (char) getc(file);
967  }
968 
969 
970  int r;
971  r = fscanf(file,"%d\n",&n);
972  if(r != 1) hit_errInternal(__func__, "Error reading CSR format", "", __FILE__, __LINE__ );
973 
974 
975  hit_cShapeNvertices(shape) = n;
976  // @arturo Ago 2015: New allocP interface
977  // hit_malloc(hit_cShapeXadj(shape), sizeof(int) * (size_t) (n+1),int*);
978  hit_malloc(hit_cShapeXadj(shape), int, n+1 );
979 
980  int i;
981  for(i=0;i<n+1;i++){
982  r = fscanf(file,"%d",& (hit_cShapeXadj(shape)[i]) );
983  if(r != 1) hit_errInternal(__func__, "Error reading CSR format", "", __FILE__, __LINE__ );
984  }
985 
986  int nedges = hit_cShapeXadj(shape)[n];
987  // @arturo Ago 2015: New allocP interface
988  // hit_malloc(hit_cShapeAdjncy(shape), sizeof(int) * (size_t) nedges, int*);
989  hit_malloc(hit_cShapeAdjncy(shape), int, nedges );
990 
991  for(i=0;i<nedges;i++){
992  r = fscanf(file,"%d",& (hit_cShapeAdjncy(shape)[i]) );
993  if(r != 1) hit_errInternal(__func__, "Error reading CSR format", "", __FILE__, __LINE__ );
994  }
995 
997  // second list
998 
999 
1000  fclose(file);
1001 
1002  return shape;
1003 }
1004 
1005 
1006 
1007 
1008 void hit_fileCSRWriteInternal(const char * csrfile, HitShape shape, int rank, const char * cfile, int line){
1009 
1010  // If the current processor is not root return a null shape.
1011  if(rank != 0){
1012  hit_warnInternal(__FUNCTION__,"Only root processor should use this function","",cfile,line);
1013  return;
1014  }
1015 
1016  FILE * file = fopen(csrfile,"w");
1017  if(file == NULL) return;
1018 
1019 
1020  fprintf(file,"# CSR graph generated by the Hitmap library\n");
1021  fprintf(file,"# Format:\n");
1022  fprintf(file,"# - Number of vertices.\n");
1023  fprintf(file,"# - Xadj list.\n");
1024  fprintf(file,"# - Adjncy list.\n");
1025  fprintf(file,"# - Vertices name list.\n\n");
1026 
1027 
1028  // Print the number of vertices
1029  int n = hit_cShapeNvertices(shape);
1030  fprintf(file,"%d\n",n);
1031 
1032  // Print xadj.
1033  int i;
1034  for(i=0;i<n+1;i++){
1035  fprintf(file,"%d ",hit_cShapeXadj(shape)[i]);
1036  }
1037  fprintf(file,"\n");
1038 
1039  // Print adjncy
1040  int nedges = hit_cShapeXadj(shape)[n];
1041  for(i=0;i<nedges;i++){
1042  fprintf(file,"%d ",hit_cShapeAdjncy(shape)[i]);
1043  }
1044  fprintf(file,"\n");
1045 
1046 
1047  fclose(file);
1048 
1049 }
#define hit_bShapeSet(bitshape, i, j)
#define hit_cShapeAdjncy(shape)
void hit_csr_renameIndexes(int Nr, int *xadj, int *adjncy)
Definition: hit_file.c:72
#define hit_tileNewType(baseType)
Definition: hit_tile.h:163
#define hit_cShapeNedges(shape)
#define hit_tileFill(var, value)
Definition: hit_tile.h:390
int hit_csr_renameIndexes_and_Symmetrize(int *pNr, int *pNc, int nz, int **pxadj, int **padjncy, const char *file, int line)
Definition: hit_file.c:86
#define HIT_ERROR
Definition: hit_error.h:54
#define hit_bShapeEdgeTarget(s, edge)
#define hit_bShapeEdgeIterator(var, shape, vertex)
Definition: hit_bshape.h:494
#define hit_calloc(ptr, type, nmemb)
Definition: hit_allocP.h:114
HitShape hit_bitmapShape(int nvertices)
Definition: hit_bshape.c:46
#define hit_cShapeCard(shape, dim)
Definition: hit_cshape.h:108
void hit_fileHBWriteBitmapInternal(const char *hbfile, HitShape shape, int rank, const char *file, int line)
Definition: hit_file.c:664
int rank
Definition: SWpar_ref.c:181
#define r(a, b, c)
int idxtype
Definition: struct.h:19
#define HIT_BITMAP_SHAPE
Definition: hit_shape.h:183
void hit_nameListCreate(HitNameList *list, int nelems)
Definition: hit_shape.c:301
#define hit_bShapeVertexIterator(var, shape)
Definition: hit_bshape.h:455
#define hit_bShapeNedges(shape)
#define hit_cShapeXadj(shape)
int hit_fileHBReadDenseInternal(const char *hbfile, int rank, void *tileP, const char *file, int line)
Definition: hit_file.c:816
#define hit_cShapeNvertices(shape)
Hitmap functions to allocate memory.
int readHB_newmat_double(const char *filename, int *M, int *N, int *nonzeros, int **colptr, int **rowind, double **val)
Definition: iohb.c:521
#define hit_cShapeNameList(shape, dim)
Definition: hit_cshape.h:162
#define hit_malloc(ptr, type, nmemb)
Definition: hit_allocP.h:93
#define HIT_NOT_USED(x)
Definition: hit_error.h:50
#define HIT_OK
Definition: hit_error.h:53
#define HIT_CSR_SHAPE
Definition: hit_shape.h:177
void hit_fileCSRWriteInternal(const char *csrfile, HitShape shape, int rank, const char *cfile, int line)
Definition: hit_file.c:1008
HitShape hit_bitmapShapeMatrix(int n, int m)
Definition: hit_bshape.c:69
HitShape shape
#define hit_bShapeNvertices(shape)
HitShape hit_fileHBRead_toBitmap_Internal(const char *fileName, int create_graph, int rank, const char *file, int line)
Definition: hit_file.c:466
HitShape hit_fileMMRead_toCSR_Internal(const char *fileName, int create_graph, int rank, const char *file, int line)
Definition: hit_file.c:246
void read_mm_matrix(char *fn, int *m, int *n, int *nz, int **i_idx, int **j_idx, double **a)
Definition: matrix_io.c:553
HitShape HIT_CSR_SHAPE_NULL
Definition: hit_shape.c:70
#define hit_tileIsNull(var)
Definition: hit_tile.h:241
HitShape hit_fileCSRReadInternal(const char *csrfile, int rank, const char *cfile, int line)
Definition: hit_file.c:942
#define hit_shapeType(s)
Definition: hit_shape.h:266
HitShape HIT_BITMAP_SHAPE_NULL
Definition: hit_shape.c:71
#define hit_realloc(ptr, type, nmemb)
Definition: hit_allocP.h:134
int compare_idxtype(const void *a, const void *b)
Definition: hit_file.c:66
HitShape hit_fileHBRead_toCSR_Internal(const char *fileName, int create_graph, int rank, const char *file, int line)
Definition: hit_file.c:340
int readHB_info(const char *filename, int *M, int *N, int *nz, char **Type, int *Nrhs)
Definition: iohb.c:234
#define idxmalloc
Definition: rename.h:383
int hit_fileHBVerticesInternal(const char *hbfile, int rank, const char *file, int line)
Definition: hit_file.c:778
#define hit_max(a, b)
Definition: hit_funcop.h:59
#define hit_warnInternal(routine, text, extraParam, file, numLine)
Definition: hit_error.h:69
#define hit_errInternal(routine, text, extraParam, file, numLine)
Definition: hit_error.h:63
void hit_fileHBWriteInternal(const char *hbfile, HitShape shape, int rank, const char *file, int line)
Definition: hit_file.c:573
#define NUMBERS_LINE