217 #include "Tpetra_Util_iohb.h"
226 char* substr(
const char* S,
const int pos,
const int len);
227 void upcase(
char* S);
228 void IOHBTerminate(
const char* message);
230 int readHB_info(
const char* filename,
int* M,
int* N,
int* nz,
char** Type,
250 int Ptrcrd, Indcrd, Valcrd, Rhscrd;
251 int Nrow, Ncol, Nnzero;
253 char Title[73], Key[9], Rhstype[4];
254 char Ptrfmt[17], Indfmt[17], Valfmt[21], Rhsfmt[21];
256 mat_type = (
char*)malloc(4);
257 if (mat_type == NULL) IOHBTerminate(
"Insufficient memory for mat_type\n");
259 if ((in_file = std::fopen(filename,
"r")) == NULL) {
260 std::fprintf(stderr,
"Error: Cannot open file: %s\n", filename);
264 readHB_header(in_file, Title, Key, mat_type, &Nrow, &Ncol, &Nnzero, Nrhs,
265 Ptrfmt, Indfmt, Valfmt, Rhsfmt,
266 &Ptrcrd, &Indcrd, &Valcrd, &Rhscrd, Rhstype);
267 std::fclose(in_file);
293 int readHB_header(std::FILE* in_file,
char* Title,
char* Key,
char* Type,
294 int* Nrow,
int* Ncol,
int* Nnzero,
int* Nrhs,
295 char* Ptrfmt,
char* Indfmt,
char* Valfmt,
char* Rhsfmt,
296 int* Ptrcrd,
int* Indcrd,
int* Valcrd,
int* Rhscrd,
301 int Totcrd, Neltvl, Nrhsix;
305 if (std::fgets(line, BUFSIZ, in_file) == NULL) {
306 std::fprintf(stderr,
"Error: Failed to read from file.\n");
309 if (std::sscanf(line,
"%*s") < 0)
310 IOHBTerminate(
"Trilinos_Util_iohb.cpp: Null (or blank) first line of HB file.\n");
311 (void)std::sscanf(line,
"%72c%8[^\n]", Title, Key);
313 *(Title + 72) =
'\0';
316 if (std::fgets(line, BUFSIZ, in_file) == NULL) {
317 std::fprintf(stderr,
"Error: Failed to read from file.\n");
320 if (std::sscanf(line,
"%*s") < 0)
321 IOHBTerminate(
"Trilinos_Util_iohb.cpp: Null (or blank) second line of HB file.\n");
322 if (std::sscanf(line,
"%i", &Totcrd) != 1) Totcrd = 0;
323 if (std::sscanf(line,
"%*i%i", Ptrcrd) != 1) *Ptrcrd = 0;
324 if (std::sscanf(line,
"%*i%*i%i", Indcrd) != 1) *Indcrd = 0;
325 if (std::sscanf(line,
"%*i%*i%*i%i", Valcrd) != 1) *Valcrd = 0;
326 if (std::sscanf(line,
"%*i%*i%*i%*i%i", Rhscrd) != 1) *Rhscrd = 0;
329 if (std::fgets(line, BUFSIZ, in_file) == NULL) {
330 std::fprintf(stderr,
"Error: Failed to read from file.\n");
333 if (std::sscanf(line,
"%*s") < 0)
334 IOHBTerminate(
"Trilinos_Util_iohb.cpp: Null (or blank) third line of HB file.\n");
335 if (std::sscanf(line,
"%3c", Type) != 1)
336 IOHBTerminate(
"Trilinos_Util_iohb.cpp: Invalid Type info, line 3 of Harwell-Boeing file.\n");
338 if (std::sscanf(line,
"%*3c%i", Nrow) != 1) *Nrow = 0;
339 if (std::sscanf(line,
"%*3c%*i%i", Ncol) != 1) *Ncol = 0;
340 if (std::sscanf(line,
"%*3c%*i%*i%i", Nnzero) != 1) *Nnzero = 0;
341 if (std::sscanf(line,
"%*3c%*i%*i%*i%i", &Neltvl) != 1) Neltvl = 0;
344 if (std::fgets(line, BUFSIZ, in_file) == NULL) {
345 std::fprintf(stderr,
"Error: Failed to read from file.\n");
348 if (std::sscanf(line,
"%*s") < 0)
349 IOHBTerminate(
"Trilinos_Util_iohb.cpp: Null (or blank) fourth line of HB file.\n");
350 if (std::sscanf(line,
"%16c", Ptrfmt) != 1)
351 IOHBTerminate(
"Trilinos_Util_iohb.cpp: Invalid format info, line 4 of Harwell-Boeing file.\n");
352 if (std::sscanf(line,
"%*16c%16c", Indfmt) != 1)
353 IOHBTerminate(
"Trilinos_Util_iohb.cpp: Invalid format info, line 4 of Harwell-Boeing file.\n");
354 if (std::sscanf(line,
"%*16c%*16c%20c", Valfmt) != 1)
355 IOHBTerminate(
"Trilinos_Util_iohb.cpp: Invalid format info, line 4 of Harwell-Boeing file.\n");
356 std::sscanf(line,
"%*16c%*16c%*20c%20c", Rhsfmt);
357 *(Ptrfmt + 16) =
'\0';
358 *(Indfmt + 16) =
'\0';
359 *(Valfmt + 20) =
'\0';
360 *(Rhsfmt + 20) =
'\0';
364 if (std::fgets(line, BUFSIZ, in_file) == NULL) {
365 std::fprintf(stderr,
"Error: Failed to read from file.\n");
368 if (std::sscanf(line,
"%*s") < 0)
369 IOHBTerminate(
"Trilinos_Util_iohb.cpp: Null (or blank) fifth line of HB file.\n");
370 if (std::sscanf(line,
"%3c", Rhstype) != 1)
371 IOHBTerminate(
"Trilinos_Util_iohb.cpp: Invalid RHS type information, line 5 of Harwell-Boeing file.\n");
372 if (std::sscanf(line,
"%*3c%i", Nrhs) != 1) *Nrhs = 0;
373 if (std::sscanf(line,
"%*3c%*i%i", &Nrhsix) != 1) Nrhsix = 0;
378 int readHB_mat_double(
const char* filename,
int colptr[],
int rowind[],
399 int i, j, ind, col, offset, count, last, Nrhs;
400 int Ptrcrd, Indcrd, Valcrd, Rhscrd;
401 int Nrow, Ncol, Nnzero, Nentries;
402 int Ptrperline, Ptrwidth, Indperline, Indwidth;
403 int Valperline, Valwidth, Valprec;
406 char Title[73], Key[9], Type[4] =
"XXX", Rhstype[4];
407 char Ptrfmt[17], Indfmt[17], Valfmt[21], Rhsfmt[21];
410 if ((in_file = std::fopen(filename,
"r")) == NULL) {
411 std::fprintf(stderr,
"Error: Cannot open file: %s\n", filename);
415 readHB_header(in_file, Title, Key, Type, &Nrow, &Ncol, &Nnzero, &Nrhs,
416 Ptrfmt, Indfmt, Valfmt, Rhsfmt,
417 &Ptrcrd, &Indcrd, &Valcrd, &Rhscrd, Rhstype);
420 ParseIfmt(Ptrfmt, &Ptrperline, &Ptrwidth);
421 ParseIfmt(Indfmt, &Indperline, &Indwidth);
422 if (Type[0] !=
'P') {
423 ParseRfmt(Valfmt, &Valperline, &Valwidth, &Valprec, &Valflag);
428 offset = 1 - _SP_base;
431 ThisElement = (
char*)malloc(Ptrwidth + 1);
432 if (ThisElement == NULL) IOHBTerminate(
"Insufficient memory for ThisElement.");
433 *(ThisElement + Ptrwidth) =
'\0';
435 for (i = 0; i < Ptrcrd; i++) {
436 if (std::fgets(line, BUFSIZ, in_file) == NULL) {
437 std::fprintf(stderr,
"Error: Failed to read from file.\n");
440 if (std::sscanf(line,
"%*s") < 0)
441 IOHBTerminate(
"Trilinos_Util_iohb.cpp: Null (or blank) line in pointer data region of HB file.\n");
443 for (ind = 0; ind < Ptrperline; ind++) {
444 if (count > Ncol)
break;
445 std::strncpy(ThisElement, line + col, Ptrwidth);
447 colptr[count] = std::atoi(ThisElement) - offset;
456 ThisElement = (
char*)malloc(Indwidth + 1);
457 if (ThisElement == NULL) IOHBTerminate(
"Insufficient memory for ThisElement.");
458 *(ThisElement + Indwidth) =
'\0';
460 for (i = 0; i < Indcrd; i++) {
461 if (std::fgets(line, BUFSIZ, in_file) == NULL) {
462 std::fprintf(stderr,
"Error: Failed to read from file.\n");
465 if (std::sscanf(line,
"%*s") < 0)
466 IOHBTerminate(
"Trilinos_Util_iohb.cpp: Null (or blank) line in index data region of HB file.\n");
468 for (ind = 0; ind < Indperline; ind++) {
469 if (count == Nnzero)
break;
470 std::strncpy(ThisElement, line + col, Indwidth);
472 rowind[count] = std::atoi(ThisElement) - offset;
481 if (Type[0] !=
'P') {
484 Nentries = 2 * Nnzero;
488 ThisElement = (
char*)malloc(Valwidth + 2);
489 if (ThisElement == NULL) IOHBTerminate(
"Insufficient memory for ThisElement.");
490 *(ThisElement + Valwidth) =
'\0';
491 *(ThisElement + Valwidth + 1) =
'\0';
493 for (i = 0; i < Valcrd; i++) {
494 if (std::fgets(line, BUFSIZ, in_file) == NULL) {
495 std::fprintf(stderr,
"Error: Failed to read from file.\n");
498 if (std::sscanf(line,
"%*s") < 0)
499 IOHBTerminate(
"Trilinos_Util_iohb.cpp: Null (or blank) line in value data region of HB file.\n");
500 if (Valflag ==
'D') {
501 while (std::strchr(line,
'D')) *std::strchr(line,
'D') =
'E';
505 for (ind = 0; ind < Valperline; ind++) {
506 if (count == Nentries)
break;
507 std::strncpy(ThisElement, line + col, Valwidth);
509 if (Valflag !=
'F' && std::strchr(ThisElement,
'E') == NULL) {
511 last = std::strlen(ThisElement);
512 for (j = last + 1; j >= 0; j--) {
513 ThisElement[j] = ThisElement[j - 1];
514 if (ThisElement[j] ==
'+' || ThisElement[j] ==
'-') {
515 ThisElement[j - 1] = Valflag;
520 val[count] = std::atof(ThisElement);
523 *(ThisElement + Valwidth) =
'\0';
524 *(ThisElement + Valwidth + 1) =
'\0';
530 std::fclose(in_file);
534 int readHB_newmat_double(
const char* filename,
int* M,
int* N,
int* nonzeros,
535 int** colptr,
int** rowind,
double** val) {
539 if (readHB_info(filename, M, N, nonzeros, &Type, &Nrhs) == 0) {
543 *colptr = (
int*)malloc((*N + 1) *
sizeof(int));
544 if (*colptr == NULL) IOHBTerminate(
"Insufficient memory for colptr.\n");
545 *rowind = (
int*)malloc(*nonzeros *
sizeof(
int));
546 if (*rowind == NULL) IOHBTerminate(
"Insufficient memory for rowind.\n");
547 if (Type[0] ==
'C') {
553 *val = (
double*)malloc(*nonzeros *
sizeof(
double) * 2);
554 if (*val == NULL) IOHBTerminate(
"Insufficient memory for val.\n");
556 if (Type[0] !=
'P') {
558 *val = (
double*)malloc(*nonzeros *
sizeof(
double));
559 if (*val == NULL) IOHBTerminate(
"Insufficient memory for val.\n");
562 return readHB_mat_double(filename, *colptr, *rowind, *val);
565 int readHB_aux_double(
const char* filename,
const char AuxType,
double b[]) {
588 int i, j, n, maxcol,
start, stride, col, last, linel;
589 int Ptrcrd, Indcrd, Valcrd, Rhscrd;
590 int Nrow, Ncol, Nnzero, Nentries;
591 int Nrhs, nvecs, rhsi;
592 int Rhsperline, Rhswidth, Rhsprec;
595 char Title[73], Key[9], Type[4] =
"XXX", Rhstype[4];
596 char Ptrfmt[17], Indfmt[17], Valfmt[21], Rhsfmt[21];
599 if ((in_file = std::fopen(filename,
"r")) == NULL) {
600 std::fprintf(stderr,
"Error: Cannot open file: %s\n", filename);
604 readHB_header(in_file, Title, Key, Type, &Nrow, &Ncol, &Nnzero, &Nrhs,
605 Ptrfmt, Indfmt, Valfmt, Rhsfmt,
606 &Ptrcrd, &Indcrd, &Valcrd, &Rhscrd, Rhstype);
609 std::fprintf(stderr,
"Warn: Attempt to read auxillary vector(s) when none are present.\n");
612 if (Rhstype[0] !=
'F') {
613 std::fprintf(stderr,
"Warn: Attempt to read auxillary vector(s) which are not stored in Full form.\n");
614 std::fprintf(stderr,
" Rhs must be specified as full. \n");
619 if (Type[0] ==
'C') {
627 if (Rhstype[1] ==
'G') nvecs++;
628 if (Rhstype[2] ==
'X') nvecs++;
630 if (AuxType ==
'G' && Rhstype[1] !=
'G') {
631 std::fprintf(stderr,
"Warn: Attempt to read auxillary Guess vector(s) when none are present.\n");
634 if (AuxType ==
'X' && Rhstype[2] !=
'X') {
635 std::fprintf(stderr,
"Warn: Attempt to read auxillary eXact solution vector(s) when none are present.\n");
639 ParseRfmt(Rhsfmt, &Rhsperline, &Rhswidth, &Rhsprec, &Rhsflag);
640 maxcol = Rhsperline * Rhswidth;
643 n = Ptrcrd + Indcrd + Valcrd;
645 for (i = 0; i < n; i++) {
646 if (std::fgets(line, BUFSIZ, in_file) == NULL) {
647 std::fprintf(stderr,
"Error: Failed to read from file.\n");
658 else if (AuxType ==
'G')
661 start = (nvecs - 1) * Nentries;
662 stride = (nvecs - 1) * Nentries;
664 if (std::fgets(line, BUFSIZ, in_file) == NULL) {
665 std::fprintf(stderr,
"Error: Failed to read from file.\n");
668 linel = std::strchr(line,
'\n') - line;
672 for (i = 0; i <
start; i++) {
673 if (col >= (maxcol < linel ? maxcol : linel)) {
674 if (std::fgets(line, BUFSIZ, in_file) == NULL) {
675 std::fprintf(stderr,
"Error: Failed to read from file.\n");
678 linel = std::strchr(line,
'\n') - line;
683 if (Rhsflag ==
'D') {
684 while (std::strchr(line,
'D')) *std::strchr(line,
'D') =
'E';
690 ThisElement = (
char*)malloc(Rhswidth + 1);
691 if (ThisElement == NULL) IOHBTerminate(
"Insufficient memory for ThisElement.");
692 *(ThisElement + Rhswidth) =
'\0';
693 for (rhsi = 0; rhsi < Nrhs; rhsi++) {
694 for (i = 0; i < Nentries; i++) {
695 if (col >= (maxcol < linel ? maxcol : linel)) {
696 if (std::fgets(line, BUFSIZ, in_file) == NULL) {
697 std::fprintf(stderr,
"Error: Failed to read from file.\n");
700 linel = std::strchr(line,
'\n') - line;
701 if (Rhsflag ==
'D') {
702 while (std::strchr(line,
'D')) *std::strchr(line,
'D') =
'E';
706 std::strncpy(ThisElement, line + col, Rhswidth);
708 if (Rhsflag !=
'F' && std::strchr(ThisElement,
'E') == NULL) {
710 last = std::strlen(ThisElement);
711 for (j = last + 1; j >= 0; j--) {
712 ThisElement[j] = ThisElement[j - 1];
713 if (ThisElement[j] ==
'+' || ThisElement[j] ==
'-') {
714 ThisElement[j - 1] = Rhsflag;
719 b[i] = std::atof(ThisElement);
725 for (i = 0; i < stride; i++) {
726 if (col >= (maxcol < linel ? maxcol : linel)) {
727 if (std::fgets(line, BUFSIZ, in_file) == NULL) {
728 std::fprintf(stderr,
"Error: Failed to read from file.\n");
731 linel = std::strchr(line,
'\n') - line;
739 std::fclose(in_file);
743 int readHB_newaux_double(
const char* filename,
const char AuxType,
double** b) {
750 readHB_info(filename, &M, &N, &nonzeros, &Type, &Nrhs);
752 std::fprintf(stderr,
"Warn: Requested read of aux vector(s) when none are present.\n");
755 if (Type[0] ==
'C') {
756 std::fprintf(stderr,
"Warning: Reading complex aux vector(s) from HB file %s.", filename);
757 std::fprintf(stderr,
" Real and imaginary parts will be interlaced in b[].");
758 *b = (
double*)malloc(M * Nrhs *
sizeof(
double) * 2);
759 if (*b == NULL) IOHBTerminate(
"Insufficient memory for rhs.\n");
760 return readHB_aux_double(filename, AuxType, *b);
762 *b = (
double*)malloc(M * Nrhs *
sizeof(
double));
763 if (*b == NULL) IOHBTerminate(
"Insufficient memory for rhs.\n");
764 return readHB_aux_double(filename, AuxType, *b);
769 int writeHB_mat_double(
const char* filename,
int M,
int N,
770 int nz,
const int colptr[],
const int rowind[],
771 const double val[],
int Nrhs,
const double rhs[],
772 const double guess[],
const double exact[],
773 const char* Title,
const char* Key,
const char* Type,
774 char* Ptrfmt,
char* Indfmt,
char* Valfmt,
char* Rhsfmt,
775 const char* Rhstype) {
786 int i, j, entry, offset, acount, linemod;
787 int totcrd, ptrcrd, indcrd, valcrd, rhscrd;
788 int nvalentries, nrhsentries;
789 int Ptrperline, Ptrwidth, Indperline, Indwidth;
790 int Rhsperline, Rhswidth, Rhsprec;
792 int Valperline, Valwidth, Valprec;
794 char pformat[16], iformat[16], vformat[19], rformat[19];
796 if (Type[0] ==
'C') {
797 nvalentries = 2 * nz;
804 if (filename != NULL) {
805 if ((out_file = std::fopen(filename,
"w")) == NULL) {
806 std::fprintf(stderr,
"Error: Cannot open file: %s\n", filename);
812 if (Ptrfmt != NULL) strcpy(Ptrfmt,
"(8I10)");
813 ParseIfmt(Ptrfmt, &Ptrperline, &Ptrwidth);
814 std::sprintf(pformat,
"%%%dd", Ptrwidth);
815 ptrcrd = (N + 1) / Ptrperline;
816 if ((N + 1) % Ptrperline != 0) ptrcrd++;
818 if (Indfmt == NULL) Indfmt = Ptrfmt;
819 ParseIfmt(Indfmt, &Indperline, &Indwidth);
820 std::sprintf(iformat,
"%%%dd", Indwidth);
821 indcrd = nz / Indperline;
822 if (nz % Indperline != 0) indcrd++;
824 if (Type[0] !=
'P') {
825 if (Valfmt != NULL) strcpy(Valfmt,
"(4E20.13)");
826 ParseRfmt(Valfmt, &Valperline, &Valwidth, &Valprec, &Valflag);
827 if (Valflag ==
'D') *std::strchr(Valfmt,
'D') =
'E';
829 std::sprintf(vformat,
"%% %d.%df", Valwidth, Valprec);
831 std::sprintf(vformat,
"%% %d.%dE", Valwidth, Valprec);
832 valcrd = nvalentries / Valperline;
833 if (nvalentries % Valperline != 0) valcrd++;
838 if (Rhsfmt == NULL) Rhsfmt = Valfmt;
839 ParseRfmt(Rhsfmt, &Rhsperline, &Rhswidth, &Rhsprec, &Rhsflag);
841 std::sprintf(rformat,
"%% %d.%df", Rhswidth, Rhsprec);
843 std::sprintf(rformat,
"%% %d.%dE", Rhswidth, Rhsprec);
844 if (Rhsflag ==
'D') *std::strchr(Rhsfmt,
'D') =
'E';
845 rhscrd = nrhsentries / Rhsperline;
846 if (nrhsentries % Rhsperline != 0) rhscrd++;
847 if (Rhstype[1] ==
'G') rhscrd += rhscrd;
848 if (Rhstype[2] ==
'X') rhscrd += rhscrd;
853 totcrd = 4 + ptrcrd + indcrd + valcrd + rhscrd;
857 std::fprintf(out_file,
"%-72s%-8s\n%14d%14d%14d%14d%14d\n", Title, Key, totcrd,
858 ptrcrd, indcrd, valcrd, rhscrd);
859 std::fprintf(out_file,
"%3s%11s%14d%14d%14d\n", Type,
" ", M, N, nz);
860 std::fprintf(out_file,
"%-16s%-16s%-20s", Ptrfmt, Indfmt, Valfmt);
864 std::fprintf(out_file,
"%-20s\n%-14s%d\n", Rhsfmt, Rhstype, Nrhs);
866 std::fprintf(out_file,
"\n");
868 offset = 1 - _SP_base;
872 for (i = 0; i < N + 1; i++) {
873 entry = colptr[i] + offset;
874 std::fprintf(out_file, pformat, entry);
875 if ((i + 1) % Ptrperline == 0) std::fprintf(out_file,
"\n");
878 if ((N + 1) % Ptrperline != 0) std::fprintf(out_file,
"\n");
881 for (i = 0; i < nz; i++) {
882 entry = rowind[i] + offset;
883 std::fprintf(out_file, iformat, entry);
884 if ((i + 1) % Indperline == 0) std::fprintf(out_file,
"\n");
887 if (nz % Indperline != 0) std::fprintf(out_file,
"\n");
891 if (Type[0] !=
'P') {
893 for (i = 0; i < nvalentries; i++) {
894 std::fprintf(out_file, vformat, val[i]);
895 if ((i + 1) % Valperline == 0) std::fprintf(out_file,
"\n");
898 if (nvalentries % Valperline != 0) std::fprintf(out_file,
"\n");
905 for (i = 0; i < Nrhs; i++) {
906 for (j = 0; j < nrhsentries; j++) {
907 std::fprintf(out_file, rformat, rhs[j]);
908 if (acount++ % Rhsperline == linemod) std::fprintf(out_file,
"\n");
910 if ((acount - 1) % Rhsperline != linemod) {
911 std::fprintf(out_file,
"\n");
912 linemod = (acount - 1) % Rhsperline;
915 if (Rhstype[1] ==
'G') {
916 for (j = 0; j < nrhsentries; j++) {
917 std::fprintf(out_file, rformat, guess[j]);
918 if (acount++ % Rhsperline == linemod) std::fprintf(out_file,
"\n");
920 if ((acount - 1) % Rhsperline != linemod) {
921 std::fprintf(out_file,
"\n");
922 linemod = (acount - 1) % Rhsperline;
924 guess += nrhsentries;
926 if (Rhstype[2] ==
'X') {
927 for (j = 0; j < nrhsentries; j++) {
928 std::fprintf(out_file, rformat, exact[j]);
929 if (acount++ % Rhsperline == linemod) std::fprintf(out_file,
"\n");
931 if ((acount - 1) % Rhsperline != linemod) {
932 std::fprintf(out_file,
"\n");
933 linemod = (acount - 1) % Rhsperline;
935 exact += nrhsentries;
941 if (std::fclose(out_file) != 0) {
942 std::fprintf(stderr,
"Error closing file in writeHB_mat_double().\n");
948 int readHB_mat_char(
const char* filename,
int colptr[],
int rowind[],
949 char val[],
char* Valfmt) {
969 int i, j, ind, col, offset, count, last;
970 int Nrow, Ncol, Nnzero, Nentries, Nrhs;
971 int Ptrcrd, Indcrd, Valcrd, Rhscrd;
972 int Ptrperline, Ptrwidth, Indperline, Indwidth;
973 int Valperline, Valwidth, Valprec;
977 char Title[73], Key[9], Type[4] =
"XXX", Rhstype[4];
978 char Ptrfmt[17], Indfmt[17], Rhsfmt[21];
980 if ((in_file = std::fopen(filename,
"r")) == NULL) {
981 std::fprintf(stderr,
"Error: Cannot open file: %s\n", filename);
985 readHB_header(in_file, Title, Key, Type, &Nrow, &Ncol, &Nnzero, &Nrhs,
986 Ptrfmt, Indfmt, Valfmt, Rhsfmt,
987 &Ptrcrd, &Indcrd, &Valcrd, &Rhscrd, Rhstype);
990 ParseIfmt(Ptrfmt, &Ptrperline, &Ptrwidth);
991 ParseIfmt(Indfmt, &Indperline, &Indwidth);
992 if (Type[0] !=
'P') {
993 ParseRfmt(Valfmt, &Valperline, &Valwidth, &Valprec, &Valflag);
994 if (Valflag ==
'D') {
995 *std::strchr(Valfmt,
'D') =
'E';
1001 offset = 1 - _SP_base;
1004 ThisElement = (
char*)malloc(Ptrwidth + 1);
1005 if (ThisElement == NULL) IOHBTerminate(
"Insufficient memory for ThisElement.");
1006 *(ThisElement + Ptrwidth) =
'\0';
1008 for (i = 0; i < Ptrcrd; i++) {
1009 if (std::fgets(line, BUFSIZ, in_file) == NULL) {
1010 std::fprintf(stderr,
"Error: Failed to read from file.\n");
1013 if (std::sscanf(line,
"%*s") < 0)
1014 IOHBTerminate(
"Trilinos_Util_iohb.cpp: Null (or blank) line in pointer data region of HB file.\n");
1016 for (ind = 0; ind < Ptrperline; ind++) {
1017 if (count > Ncol)
break;
1018 std::strncpy(ThisElement, line + col, Ptrwidth);
1020 colptr[count] = std::atoi(ThisElement) - offset;
1029 ThisElement = (
char*)malloc(Indwidth + 1);
1030 if (ThisElement == NULL) IOHBTerminate(
"Insufficient memory for ThisElement.");
1031 *(ThisElement + Indwidth) =
'\0';
1033 for (i = 0; i < Indcrd; i++) {
1034 if (std::fgets(line, BUFSIZ, in_file) == NULL) {
1035 std::fprintf(stderr,
"Error: Failed to read from file.\n");
1038 if (std::sscanf(line,
"%*s") < 0)
1039 IOHBTerminate(
"Trilinos_Util_iohb.cpp: Null (or blank) line in index data region of HB file.\n");
1041 for (ind = 0; ind < Indperline; ind++) {
1042 if (count == Nnzero)
break;
1043 std::strncpy(ThisElement, line + col, Indwidth);
1045 rowind[count] = std::atoi(ThisElement) - offset;
1054 if (Type[0] !=
'P') {
1057 Nentries = 2 * Nnzero;
1061 ThisElement = (
char*)malloc(Valwidth + 1);
1062 if (ThisElement == NULL) IOHBTerminate(
"Insufficient memory for ThisElement.");
1063 *(ThisElement + Valwidth) =
'\0';
1065 for (i = 0; i < Valcrd; i++) {
1066 if (std::fgets(line, BUFSIZ, in_file) == NULL) {
1067 std::fprintf(stderr,
"Error: Failed to read from file.\n");
1070 if (std::sscanf(line,
"%*s") < 0)
1071 IOHBTerminate(
"Trilinos_Util_iohb.cpp: Null (or blank) line in value data region of HB file.\n");
1072 if (Valflag ==
'D') {
1073 while (std::strchr(line,
'D')) *std::strchr(line,
'D') =
'E';
1076 for (ind = 0; ind < Valperline; ind++) {
1077 if (count == Nentries)
break;
1078 ThisElement = &val[count * Valwidth];
1079 std::strncpy(ThisElement, line + col, Valwidth);
1081 if (Valflag !=
'F' && std::strchr(ThisElement,
'E') == NULL) {
1083 last = std::strlen(ThisElement);
1084 for (j = last + 1; j >= 0; j--) {
1085 ThisElement[j] = ThisElement[j - 1];
1086 if (ThisElement[j] ==
'+' || ThisElement[j] ==
'-') {
1087 ThisElement[j - 1] = Valflag;
1101 int readHB_newmat_char(
const char* filename,
int* M,
int* N,
int* nonzeros,
int** colptr,
1102 int** rowind,
char** val,
char** Valfmt) {
1105 int Ptrcrd, Indcrd, Valcrd, Rhscrd;
1106 int Valperline, Valwidth, Valprec;
1108 char Title[73], Key[9], Type[4] =
"XXX", Rhstype[4];
1109 char Ptrfmt[17], Indfmt[17], Rhsfmt[21];
1111 if ((in_file = std::fopen(filename,
"r")) == NULL) {
1112 std::fprintf(stderr,
"Error: Cannot open file: %s\n", filename);
1116 *Valfmt = (
char*)malloc(21 *
sizeof(
char));
1117 if (*Valfmt == NULL) IOHBTerminate(
"Insufficient memory for Valfmt.");
1118 readHB_header(in_file, Title, Key, Type, M, N, nonzeros, &Nrhs,
1119 Ptrfmt, Indfmt, (*Valfmt), Rhsfmt,
1120 &Ptrcrd, &Indcrd, &Valcrd, &Rhscrd, Rhstype);
1121 std::fclose(in_file);
1122 ParseRfmt(*Valfmt, &Valperline, &Valwidth, &Valprec, &Valflag);
1124 *colptr = (
int*)malloc((*N + 1) *
sizeof(int));
1125 if (*colptr == NULL) IOHBTerminate(
"Insufficient memory for colptr.\n");
1126 *rowind = (
int*)malloc(*nonzeros *
sizeof(
int));
1127 if (*rowind == NULL) IOHBTerminate(
"Insufficient memory for rowind.\n");
1128 if (Type[0] ==
'C') {
1134 *val = (
char*)malloc(*nonzeros * Valwidth *
sizeof(
char) * 2);
1135 if (*val == NULL) IOHBTerminate(
"Insufficient memory for val.\n");
1137 if (Type[0] !=
'P') {
1139 *val = (
char*)malloc(*nonzeros * Valwidth *
sizeof(
char));
1140 if (*val == NULL) IOHBTerminate(
"Insufficient memory for val.\n");
1143 return readHB_mat_char(filename, *colptr, *rowind, *val, *Valfmt);
1146 int readHB_aux_char(
const char* filename,
const char AuxType,
char b[]) {
1169 int i, j, n, maxcol,
start, stride, col, last, linel, nvecs, rhsi;
1170 int Nrow, Ncol, Nnzero, Nentries, Nrhs;
1171 int Ptrcrd, Indcrd, Valcrd, Rhscrd;
1172 int Rhsperline, Rhswidth, Rhsprec;
1174 char Title[73], Key[9], Type[4] =
"XXX", Rhstype[4];
1175 char Ptrfmt[17], Indfmt[17], Valfmt[21], Rhsfmt[21];
1179 if ((in_file = std::fopen(filename,
"r")) == NULL) {
1180 std::fprintf(stderr,
"Error: Cannot open file: %s\n", filename);
1184 readHB_header(in_file, Title, Key, Type, &Nrow, &Ncol, &Nnzero, &Nrhs,
1185 Ptrfmt, Indfmt, Valfmt, Rhsfmt,
1186 &Ptrcrd, &Indcrd, &Valcrd, &Rhscrd, Rhstype);
1189 std::fprintf(stderr,
"Warn: Attempt to read auxillary vector(s) when none are present.\n");
1192 if (Rhstype[0] !=
'F') {
1193 std::fprintf(stderr,
"Warn: Attempt to read auxillary vector(s) which are not stored in Full form.\n");
1194 std::fprintf(stderr,
" Rhs must be specified as full. \n");
1199 if (Type[0] ==
'C') {
1200 Nentries = 2 * Nrow;
1207 if (Rhstype[1] ==
'G') nvecs++;
1208 if (Rhstype[2] ==
'X') nvecs++;
1210 if (AuxType ==
'G' && Rhstype[1] !=
'G') {
1211 std::fprintf(stderr,
"Warn: Attempt to read auxillary Guess vector(s) when none are present.\n");
1214 if (AuxType ==
'X' && Rhstype[2] !=
'X') {
1215 std::fprintf(stderr,
"Warn: Attempt to read auxillary eXact solution vector(s) when none are present.\n");
1219 ParseRfmt(Rhsfmt, &Rhsperline, &Rhswidth, &Rhsprec, &Rhsflag);
1220 maxcol = Rhsperline * Rhswidth;
1223 n = Ptrcrd + Indcrd + Valcrd;
1225 for (i = 0; i < n; i++) {
1226 if (std::fgets(line, BUFSIZ, in_file) == NULL) {
1227 std::fprintf(stderr,
"Error: Failed to read from file.\n");
1238 else if (AuxType ==
'G')
1241 start = (nvecs - 1) * Nentries;
1242 stride = (nvecs - 1) * Nentries;
1244 if (std::fgets(line, BUFSIZ, in_file) == NULL) {
1245 std::fprintf(stderr,
"Error: Failed to read from file.\n");
1248 linel = std::strchr(line,
'\n') - line;
1249 if (std::sscanf(line,
"%*s") < 0)
1250 IOHBTerminate(
"Trilinos_Util_iohb.cpp: Null (or blank) line in auxillary vector data region of HB file.\n");
1254 for (i = 0; i <
start; i++) {
1256 if (col >= (maxcol < linel ? maxcol : linel)) {
1257 if (std::fgets(line, BUFSIZ, in_file) == NULL) {
1258 std::fprintf(stderr,
"Error: Failed to read from file.\n");
1261 linel = std::strchr(line,
'\n') - line;
1262 if (std::sscanf(line,
"%*s") < 0)
1263 IOHBTerminate(
"Trilinos_Util_iohb.cpp: Null (or blank) line in auxillary vector data region of HB file.\n");
1268 if (Rhsflag ==
'D') {
1269 while (std::strchr(line,
'D')) *std::strchr(line,
'D') =
'E';
1274 for (rhsi = 0; rhsi < Nrhs; rhsi++) {
1275 for (i = 0; i < Nentries; i++) {
1276 if (col >= (maxcol < linel ? maxcol : linel)) {
1277 if (std::fgets(line, BUFSIZ, in_file) == NULL) {
1278 std::fprintf(stderr,
"Error: Failed to read from file.\n");
1281 linel = std::strchr(line,
'\n') - line;
1282 if (std::sscanf(line,
"%*s") < 0)
1283 IOHBTerminate(
"Trilinos_Util_iohb.cpp: Null (or blank) line in auxillary vector data region of HB file.\n");
1284 if (Rhsflag ==
'D') {
1285 while (std::strchr(line,
'D')) *std::strchr(line,
'D') =
'E';
1289 ThisElement = &b[i * Rhswidth];
1290 std::strncpy(ThisElement, line + col, Rhswidth);
1291 if (Rhsflag !=
'F' && std::strchr(ThisElement,
'E') == NULL) {
1293 last = std::strlen(ThisElement);
1294 for (j = last + 1; j >= 0; j--) {
1295 ThisElement[j] = ThisElement[j - 1];
1296 if (ThisElement[j] ==
'+' || ThisElement[j] ==
'-') {
1297 ThisElement[j - 1] = Rhsflag;
1304 b += Nentries * Rhswidth;
1308 for (i = 0; i < stride; i++) {
1310 if (col >= (maxcol < linel ? maxcol : linel)) {
1311 if (std::fgets(line, BUFSIZ, in_file) == NULL) {
1312 std::fprintf(stderr,
"Error: Failed to read from file.\n");
1315 linel = std::strchr(line,
'\n') - line;
1316 if (std::sscanf(line,
"%*s") < 0)
1317 IOHBTerminate(
"Trilinos_Util_iohb.cpp: Null (or blank) line in auxillary vector data region of HB file.\n");
1323 std::fclose(in_file);
1327 int readHB_newaux_char(
const char* filename,
const char AuxType,
char** b,
char** Rhsfmt) {
1329 int Ptrcrd, Indcrd, Valcrd, Rhscrd;
1330 int Nrow, Ncol, Nnzero, Nrhs;
1331 int Rhsperline, Rhswidth, Rhsprec;
1333 char Title[73], Key[9], Type[4] =
"XXX", Rhstype[4];
1334 char Ptrfmt[17], Indfmt[17], Valfmt[21];
1336 if ((in_file = std::fopen(filename,
"r")) == NULL) {
1337 std::fprintf(stderr,
"Error: Cannot open file: %s\n", filename);
1341 *Rhsfmt = (
char*)malloc(21 *
sizeof(
char));
1342 if (*Rhsfmt == NULL) IOHBTerminate(
"Insufficient memory for Rhsfmt.");
1343 readHB_header(in_file, Title, Key, Type, &Nrow, &Ncol, &Nnzero, &Nrhs,
1344 Ptrfmt, Indfmt, Valfmt, (*Rhsfmt),
1345 &Ptrcrd, &Indcrd, &Valcrd, &Rhscrd, Rhstype);
1346 std::fclose(in_file);
1348 std::fprintf(stderr,
"Warn: Requested read of aux vector(s) when none are present.\n");
1351 ParseRfmt(*Rhsfmt, &Rhsperline, &Rhswidth, &Rhsprec, &Rhsflag);
1352 if (Type[0] ==
'C') {
1353 std::fprintf(stderr,
"Warning: Reading complex aux vector(s) from HB file %s.", filename);
1354 std::fprintf(stderr,
" Real and imaginary parts will be interlaced in b[].");
1355 *b = (
char*)malloc(Nrow * Nrhs * Rhswidth *
sizeof(
char) * 2);
1356 if (*b == NULL) IOHBTerminate(
"Insufficient memory for rhs.\n");
1357 return readHB_aux_char(filename, AuxType, *b);
1359 *b = (
char*)malloc(Nrow * Nrhs * Rhswidth *
sizeof(
char));
1360 if (*b == NULL) IOHBTerminate(
"Insufficient memory for rhs.\n");
1361 return readHB_aux_char(filename, AuxType, *b);
1366 int writeHB_mat_char(
const char* filename,
int M,
int N,
1367 int nz,
const int colptr[],
const int rowind[],
1368 const char val[],
int Nrhs,
const char rhs[],
1369 const char guess[],
const char exact[],
1370 const char* Title,
const char* Key,
const char* Type,
1371 char* Ptrfmt,
char* Indfmt,
char* Valfmt,
char* Rhsfmt,
1372 const char* Rhstype) {
1382 std::FILE* out_file;
1383 int i, j, acount, linemod, entry, offset;
1384 int totcrd, ptrcrd, indcrd, valcrd, rhscrd;
1385 int nvalentries, nrhsentries;
1386 int Ptrperline, Ptrwidth, Indperline, Indwidth;
1387 int Rhsperline, Rhswidth, Rhsprec;
1389 int Valperline, Valwidth, Valprec;
1391 char pformat[16], iformat[16], vformat[19], rformat[19];
1393 if (Type[0] ==
'C') {
1394 nvalentries = 2 * nz;
1395 nrhsentries = 2 * M;
1401 if (filename != NULL) {
1402 if ((out_file = std::fopen(filename,
"w")) == NULL) {
1403 std::fprintf(stderr,
"Error: Cannot open file: %s\n", filename);
1409 if (Ptrfmt != NULL) strcpy(Ptrfmt,
"(8I10)");
1410 ParseIfmt(Ptrfmt, &Ptrperline, &Ptrwidth);
1411 std::sprintf(pformat,
"%%%dd", Ptrwidth);
1413 if (Indfmt == NULL) Indfmt = Ptrfmt;
1414 ParseIfmt(Indfmt, &Indperline, &Indwidth);
1415 std::sprintf(iformat,
"%%%dd", Indwidth);
1417 if (Type[0] !=
'P') {
1418 if (Valfmt != NULL) strcpy(Valfmt,
"(4E20.13)");
1419 ParseRfmt(Valfmt, &Valperline, &Valwidth, &Valprec, &Valflag);
1420 std::sprintf(vformat,
"%%%ds", Valwidth);
1423 ptrcrd = (N + 1) / Ptrperline;
1424 if ((N + 1) % Ptrperline != 0) ptrcrd++;
1426 indcrd = nz / Indperline;
1427 if (nz % Indperline != 0) indcrd++;
1429 valcrd = nvalentries / Valperline;
1430 if (nvalentries % Valperline != 0) valcrd++;
1433 if (Rhsfmt == NULL) Rhsfmt = Valfmt;
1434 ParseRfmt(Rhsfmt, &Rhsperline, &Rhswidth, &Rhsprec, &Rhsflag);
1435 std::sprintf(rformat,
"%%%ds", Rhswidth);
1436 rhscrd = nrhsentries / Rhsperline;
1437 if (nrhsentries % Rhsperline != 0) rhscrd++;
1438 if (Rhstype[1] ==
'G') rhscrd += rhscrd;
1439 if (Rhstype[2] ==
'X') rhscrd += rhscrd;
1444 totcrd = 4 + ptrcrd + indcrd + valcrd + rhscrd;
1448 std::fprintf(out_file,
"%-72s%-8s\n%14d%14d%14d%14d%14d\n", Title, Key, totcrd,
1449 ptrcrd, indcrd, valcrd, rhscrd);
1450 std::fprintf(out_file,
"%3s%11s%14d%14d%14d\n", Type,
" ", M, N, nz);
1451 std::fprintf(out_file,
"%-16s%-16s%-20s", Ptrfmt, Indfmt, Valfmt);
1455 std::fprintf(out_file,
"%-20s\n%-14s%d\n", Rhsfmt, Rhstype, Nrhs);
1457 std::fprintf(out_file,
"\n");
1459 offset = 1 - _SP_base;
1463 for (i = 0; i < N + 1; i++) {
1464 entry = colptr[i] + offset;
1465 std::fprintf(out_file, pformat, entry);
1466 if ((i + 1) % Ptrperline == 0) std::fprintf(out_file,
"\n");
1469 if ((N + 1) % Ptrperline != 0) std::fprintf(out_file,
"\n");
1472 for (i = 0; i < nz; i++) {
1473 entry = rowind[i] + offset;
1474 std::fprintf(out_file, iformat, entry);
1475 if ((i + 1) % Indperline == 0) std::fprintf(out_file,
"\n");
1478 if (nz % Indperline != 0) std::fprintf(out_file,
"\n");
1482 if (Type[0] !=
'P') {
1483 for (i = 0; i < nvalentries; i++) {
1484 std::fprintf(out_file, vformat, val + i * Valwidth);
1485 if ((i + 1) % Valperline == 0) std::fprintf(out_file,
"\n");
1488 if (nvalentries % Valperline != 0) std::fprintf(out_file,
"\n");
1494 for (j = 0; j < Nrhs; j++) {
1495 for (i = 0; i < nrhsentries; i++) {
1496 std::fprintf(out_file, rformat, rhs + i * Rhswidth);
1497 if (acount++ % Rhsperline == linemod) std::fprintf(out_file,
"\n");
1499 if (acount % Rhsperline != linemod) {
1500 std::fprintf(out_file,
"\n");
1501 linemod = (acount - 1) % Rhsperline;
1503 if (Rhstype[1] ==
'G') {
1504 for (i = 0; i < nrhsentries; i++) {
1505 std::fprintf(out_file, rformat, guess + i * Rhswidth);
1506 if (acount++ % Rhsperline == linemod) std::fprintf(out_file,
"\n");
1508 if (acount % Rhsperline != linemod) {
1509 std::fprintf(out_file,
"\n");
1510 linemod = (acount - 1) % Rhsperline;
1513 if (Rhstype[2] ==
'X') {
1514 for (i = 0; i < nrhsentries; i++) {
1515 std::fprintf(out_file, rformat, exact + i * Rhswidth);
1516 if (acount++ % Rhsperline == linemod) std::fprintf(out_file,
"\n");
1518 if (acount % Rhsperline != linemod) {
1519 std::fprintf(out_file,
"\n");
1520 linemod = (acount - 1) % Rhsperline;
1527 if (std::fclose(out_file) != 0) {
1528 std::fprintf(stderr,
"Error closing file in writeHB_mat_char().\n");
1534 int ParseIfmt(
char* fmt,
int* perline,
int* width) {
1546 tmp = std::strchr(fmt,
'(');
1547 tmp = substr(fmt, tmp - fmt + 1, std::strchr(fmt,
'I') - tmp - 1);
1548 *perline = std::atoi(tmp);
1549 if (*perline == 0) *perline = 1;
1550 if (tmp != NULL) free((
void*)tmp);
1551 tmp = std::strchr(fmt,
'I');
1552 tmp = substr(fmt, tmp - fmt + 1, std::strchr(fmt,
')') - tmp - 1);
1553 *width = std::atoi(tmp);
1554 if (tmp != NULL) free((
void*)tmp);
1558 int ParseRfmt(
char* fmt,
int* perline,
int* width,
int* prec,
int* flag) {
1579 if (std::strchr(fmt,
'(') != NULL) fmt = std::strchr(fmt,
'(');
1580 if (std::strchr(fmt,
')') != NULL) {
1581 tmp2 = std::strchr(fmt,
')');
1582 while (std::strchr(tmp2 + 1,
')') != NULL) {
1583 tmp2 = std::strchr(tmp2 + 1,
')');
1587 if (std::strchr(fmt,
'P') != NULL)
1589 if (std::strchr(fmt,
'(') != NULL) {
1590 tmp = std::strchr(fmt,
'P');
1591 if (*(++tmp) ==
',') tmp++;
1592 tmp3 = std::strchr(fmt,
'(') + 1;
1595 while (*(tmp2 + len) !=
'\0') {
1596 *tmp2 = *(tmp2 + len);
1599 *(std::strchr(fmt,
')') + 1) =
'\0';
1602 if (std::strchr(fmt,
'E') != NULL) {
1604 }
else if (std::strchr(fmt,
'D') != NULL) {
1606 }
else if (std::strchr(fmt,
'F') != NULL) {
1609 std::fprintf(stderr,
"Real format %s in H/B file not supported.\n", fmt);
1612 tmp = std::strchr(fmt,
'(');
1613 tmp = substr(fmt, tmp - fmt + 1, std::strchr(fmt, *flag) - tmp - 1);
1614 *perline = std::atoi(tmp);
1615 if (*perline == 0) *perline = 1;
1616 if (tmp != NULL) free((
void*)tmp);
1617 tmp = std::strchr(fmt, *flag);
1618 if (std::strchr(fmt,
'.')) {
1619 tmp1 = substr(fmt, std::strchr(fmt,
'.') - fmt + 1, std::strchr(fmt,
')') - std::strchr(fmt,
'.') - 1);
1620 *prec = std::atoi(tmp1);
1621 if (tmp1 != NULL) free((
void*)tmp1);
1622 tmp1 = substr(fmt, tmp - fmt + 1, std::strchr(fmt,
'.') - tmp - 1);
1624 tmp1 = substr(fmt, tmp - fmt + 1, std::strchr(fmt,
')') - tmp - 1);
1626 *width = std::atoi(tmp1);
1627 if (tmp1 != NULL) free((
void*)tmp1);
1631 char* substr(
const char* S,
const int pos,
const int len) {
1634 if ((
size_t)pos + len <= std::strlen(S)) {
1635 SubS = (
char*)malloc(len + 1);
1636 if (SubS == NULL) IOHBTerminate(
"Insufficient memory for SubS.");
1637 for (i = 0; i < len; i++) SubS[i] = S[pos + i];
1646 void upcase(
char* S) {
1649 len = std::strlen(S);
1650 for (i = 0; i < len; i++)
1651 S[i] = std::toupper(S[i]);
1654 void IOHBTerminate(
const char* message) {
1655 std::fprintf(stderr,
"%s", message);
void start()
Start the deep_copy counter.