1: #include "mmloader.h"

  3: PetscErrorCode MatCreateFromMTX(Mat *A, const char *filein, PetscBool aijonly)
  4: {
  5:   MM_typecode  matcode;
  6:   FILE        *file;
  7:   PetscInt     M, N, ninput;
  8:   PetscInt    *ia, *ja;
  9:   PetscInt     i, j, nz, *rownz;
 10:   PetscScalar *val;
 11:   PetscBool    sametype, symmetric = PETSC_FALSE, skew = PETSC_FALSE;

 13:   /* Read in matrix */
 14:   PetscFunctionBeginUser;
 15:   PetscCall(PetscFOpen(PETSC_COMM_SELF, filein, "r", &file));
 16:   PetscCheck(mm_read_banner(file, &matcode) == 0, PETSC_COMM_SELF, PETSC_ERR_ARG_WRONG, "Could not process Matrix Market banner.");
 17:   /*  This is how one can screen matrix types if their application */
 18:   /*  only supports a subset of the Matrix Market data types.      */
 19:   PetscCheck(mm_is_matrix(matcode) && mm_is_sparse(matcode) && (mm_is_real(matcode) || mm_is_integer(matcode)), PETSC_COMM_SELF, PETSC_ERR_ARG_WRONG, "Input must be a sparse real or integer matrix. Market Market type: [%s]", mm_typecode_to_str(matcode));

 21:   if (mm_is_symmetric(matcode)) symmetric = PETSC_TRUE;
 22:   if (mm_is_skew(matcode)) skew = PETSC_TRUE;

 24:   /* Find out size of sparse matrix .... mmio.h declares these as plain int, so
 25:      read into int temporaries: PetscInt* != int* under 64-bit indices. */
 26:   {
 27:     int Mi, Ni, nzi;

 29:     PetscCheck(mm_read_mtx_crd_size(file, &Mi, &Ni, &nzi) == 0, PETSC_COMM_SELF, PETSC_ERR_ARG_WRONG, "Size of sparse matrix is wrong.");
 30:     M  = Mi;
 31:     N  = Ni;
 32:     nz = nzi;
 33:   }

 35:   /* Reserve memory for matrices */
 36:   PetscCall(PetscMalloc4(nz, &ia, nz, &ja, nz, &val, M, &rownz));
 37:   for (i = 0; i < M; i++) rownz[i] = 0;

 39:   /* NOTE: when reading in doubles, ANSI C requires the use of the "l"  */
 40:   /*   specifier as in "%lg", "%lf", "%le", otherwise errors will occur */
 41:   /*  (ANSI C X3.159-1989, Sec. 4.9.6.2, p. 136 lines 13-15)            */
 42:   for (i = 0; i < nz; i++) {
 43:     int row, col; /* %d reads int; ia[]/ja[] are PetscInt (long) under 64-bit indices */

 45:     ninput = fscanf(file, "%d %d %lg\n", &row, &col, &val[i]);
 46:     PetscCheck(ninput >= 3, PETSC_COMM_SELF, PETSC_ERR_FILE_UNEXPECTED, "Badly formatted input file");
 47:     ia[i] = row - 1;
 48:     ja[i] = col - 1;                      /* adjust from 1-based to 0-based */
 49:     if ((symmetric && aijonly) || skew) { /* transpose */
 50:       rownz[ia[i]]++;
 51:       if (ja[i] != ia[i]) rownz[ja[i]]++;
 52:     } else {
 53:       if (symmetric) rownz[ja[i]]++;
 54:       else rownz[ia[i]]++;
 55:     }
 56:   }
 57:   PetscCall(PetscFClose(PETSC_COMM_SELF, file));

 59:   /* Create, preallocate, and then assemble the matrix */
 60:   PetscCall(MatCreate(PETSC_COMM_SELF, A));
 61:   PetscCall(MatSetSizes(*A, PETSC_DECIDE, PETSC_DECIDE, M, N));

 63:   if (symmetric && !aijonly) {
 64:     PetscCall(MatSetType(*A, MATSEQSBAIJ));
 65:     PetscCall(MatSetFromOptions(*A));
 66:     PetscCall(MatSeqSBAIJSetPreallocation(*A, 1, 0, rownz));
 67:     PetscCall(PetscObjectTypeCompare((PetscObject)*A, MATSEQSBAIJ, &sametype));
 68:     PetscCheck(sametype, PETSC_COMM_SELF, PETSC_ERR_ARG_WRONG, "Only AIJ and SBAIJ are supported. Your mattype is not supported");
 69:   } else {
 70:     PetscCall(MatSetType(*A, MATSEQAIJ));
 71:     PetscCall(MatSetFromOptions(*A));
 72:     PetscCall(MatSeqAIJSetPreallocation(*A, 0, rownz));
 73:     PetscCall(PetscObjectTypeCompare((PetscObject)*A, MATSEQAIJ, &sametype));
 74:     PetscCheck(sametype, PETSC_COMM_SELF, PETSC_ERR_ARG_WRONG, "Only AIJ and SBAIJ are supported. Your mattype is not supported");
 75:   }
 76:   /* Add values to the matrix, these correspond to lower triangular part for symmetric or skew matrices */
 77:   if (!(symmetric && !aijonly))
 78:     for (j = 0; j < nz; j++) PetscCall(MatSetValues(*A, 1, &ia[j], 1, &ja[j], &val[j], INSERT_VALUES));

 80:   /* Add values to the upper triangular part for symmetric matrices. MatrixMarket stores a symmetric matrix in its lower triangular part, so insert its transpose. For SBAIJ this fills the stored upper triangle; for AIJ (-aij_only) it completes the full symmetric matrix. */
 81:   if (symmetric)
 82:     for (j = 0; j < nz; j++) PetscCall(MatSetValues(*A, 1, &ja[j], 1, &ia[j], &val[j], INSERT_VALUES));
 83:   if (skew) {
 84:     for (j = 0; j < nz; j++) {
 85:       val[j] = -val[j];
 86:       PetscCall(MatSetValues(*A, 1, &ja[j], 1, &ia[j], &val[j], INSERT_VALUES));
 87:     }
 88:   }
 89:   PetscCall(MatAssemblyBegin(*A, MAT_FINAL_ASSEMBLY));
 90:   PetscCall(MatAssemblyEnd(*A, MAT_FINAL_ASSEMBLY));
 91:   PetscCall(PetscFree4(ia, ja, val, rownz));
 92:   PetscFunctionReturn(PETSC_SUCCESS);
 93: }