trilinos / Trilinos

Primary repository for the Trilinos Project
https://trilinos.org/
Other
1.18k stars 559 forks source link

Amesos2: Reindexing of matrix #13205

Open maxfirmbach opened 4 days ago

maxfirmbach commented 4 days ago

Question

Dear Amesos2-developers,

currently I'm trying to switch our in-house code from Amesos to Amesos2. As, let's call it pre-processing step, we do a reindexing of our linear problem to have a continous GID numbering from 0 ... n, using EpetraExt (which in our case is necessary as we have GID jumps). Is there a similar feature implemented in Xpetra or Tpetra?

I've seen that Amesos2 has a reindex input parameter, but looking into the code, it seems to do ... basically nothing? From documentation:

We plan in the future to support the following parameters:

    "Reindex" : { true | false }. Put the matrix row and column indices into the range [0..n]. 

Are there any intentions to implement this feature properly into Amesos2?

Best regards, Max

@trilinos/amesos2 @mayrmt

csiefer2 commented 4 days ago

@ndellingwood @iyamazaki

MalachiTimothyPhillips commented 4 days ago

As an aside, I don't think you are required to have continuous GIDs so long as you set the IsContiguous setting to false. For our application code, we've never had to remap our discontinuous GIDs to continuous GIDs.

See, for example: https://github.com/trilinos/Trilinos/blob/master/packages/amesos2/test/solvers/Superlu_UnitTests.cpp#L372.

The one big caveat with that is if you dump the matrix to matrix-market format, you'll need to remap those GIDs to be 1-based index to be in compliance with the matrix market format. Here's a completely un-tested, un-supported, provided as-is without any warranty Python script that can accomplish that as a post-processing step:

Script: ```python #!/usr/bin/env python3 import os.path import scipy.io import numpy as np import sys from itertools import (takewhile, repeat) import glob import pathlib import argparse def warn_on_missing_file(fileName): errString = f"File {fileName} does not exist!" print(errString) return def warn_on_existing_file(fileName): errString = f"Cowardly refusing to overwrite file {fileName}.\n" errString += "Please use --overwrite to allow overwriting to occur, or alter the output matrix filename.\n" print(errString) return def draw_progress_bar(percentage, process="", num_bins=20): num_draw = int(percentage * num_bins) sys.stdout.write('\r') sys.stdout.write( f"{process} : [{'='*num_draw}{' ' * (num_bins-num_draw)}] {int(100 * percentage)}%") sys.stdout.flush() return def grab_line_count(filename): f = open(filename, 'rb') bufgen = takewhile(lambda x: x, (f.raw.read(1024*1024) for _ in repeat(None))) return sum(buf.count(b'\n') for buf in bufgen) def read_matrix_file(matrixFileName): gids = {} A_coo = {} readSizeLine = False update_line_freq = 1000 num_lines = grab_line_count(matrixFileName) header_content = "" with open(matrixFileName, "r") as matrixFile: for line_num, line in enumerate(matrixFile): if line_num % update_line_freq == 0: draw_progress_bar((line_num+1)/num_lines, "reading matrix file") if "%" in line: header_content += line continue if not readSizeLine: header_content += line readSizeLine = True continue i, j, A_ij = line.split() i, j = int(i), int(j) A_coo[(i, j)] = A_ij gids[i] = None draw_progress_bar(1, "reading matrix file") print() rows = {} row = 1 for gid in gids.keys(): rows[gid] = row row = row + 1 return rows, A_coo, header_content def write_matrix_file(rows, A_coo, header, outputMatrixFileName): update_line_freq = 1000 newFileContents = header num_matrix_entries = len(A_coo) for mat_entry, ((gid_i, gid_j), A_ij) in enumerate(A_coo.items()): if mat_entry % update_line_freq == 0: draw_progress_bar((mat_entry+1)/num_matrix_entries, "creating new matrix file") i = rows[gid_i] j = rows[gid_j] newFileContents += f"{i} {j} {A_ij}\n" draw_progress_bar(1, "creating new matrix file") print() with open(outputMatrixFileName, "w") as newMatrixFile: newMatrixFile.write(newFileContents) return def reassign_matrix_rows(matrixFileName, newMatrixFileName): rows, A_coo, header = read_matrix_file(matrixFileName) write_matrix_file(rows, A_coo, header, newMatrixFileName) return if __name__ == "__main__": parser = argparse.ArgumentParser( description='Re-assign row ids in matrix market file to force them to be contiguous and 1-base indexed') parser.add_argument('-i', '--input-matrix', help='Input matrix-market file', required=True) parser.add_argument('-o', '--output-matrix', help='Output matrix-market file', required=True) parser.add_argument('--overwrite', action='store_true') args = vars(parser.parse_args()) inputMatrix = args['input_matrix'] outputMatrix = args['output_matrix'] allGood = True if not os.path.isfile(inputMatrix): warn_on_missing_file(inputMatrix) print("Please check your --input-matrix/-i option.") allGood = False if os.path.isfile(outputMatrix) and not overwrite: warn_on_existing_file(outputMatrix) allGood = False reassign_matrix_rows(inputMatrix, outputMatrix) ```
maxfirmbach commented 3 days ago

@MalachiTimothyPhillips Interesting! I will try to set the IsContiguous parameter to false, maybe that already solves the problem.