Parse large ASCII matrix

I am parsing a large ASCII sparse matrix, output of PETSc, but my current implementation is impossibly slow. I don't know how I could speed that up.

The ASCII file of a sparse PETSc aij matrix is as follows (with col being an integer index, and val being a floating point, e.g., 1.372e-1):

row0: (col, val) (col, val) ...
row1: (col, val) (col, val) ...
row2: (col, val) (col, val) ...
...

What I did right now is simple: with a text editor I remove the rowX: part and then use the following code:

(* Open and split into lines *)
Ktxt = StringSplit[Import["K.txt"], "\n"]

(* convert ASCII exponents *)
convert[s_] := ToExpression /@ StringReplace[StringSplit@s, "e" | "E" :> "*^"]

l = {}

(* For each row, convert and append element *)
For[r = 1, r <= Length[Ktxt], r++,
t = StringSplit[Ktxt[[r]]];
For[i = 1, i <= Length[t], i += 2,
s1 = (convert[StringReplace[t[[i]], {"(" -> "", "," -> ""}]])[];
s2 = (convert[StringReplace[t[[i + 1]], ")" -> ""]])[];
AppendTo[l, {r, s1 + 1} -> s2]
];
If[Mod[r, 100] == 0, Print["row ", r]]
]

(* Make the sparse matrix *)
MM = SparseArray[l]

Is there a way to make this faster? For a 15,000 x 15,000 matrix with 800,000 non-zero elements it takes hours.

Thank you!

• Why not have PETSc save your matrix to a format that Mathematica can Import[], like *.mtx? May 29 '17 at 12:26
• First, I am not skilled in PETSc, which is one package used by another software I am using, and second, my goal is to improve my Mathematica knowledge. May 29 '17 at 12:43
• Please provide a sample aij matrix file (or code that generates its content) -- it is much more likely to get responses to your question. May 29 '17 at 13:06
• I modified the question: that sample code is the file itself, actually. May 29 '17 at 13:16
• cant test here, but what does Import[..,"Table"] do with that? May 29 '17 at 13:42

You are doing string operations on a much finer granularity level than necessary. The code below is fairly fast on my computer for matrices with the size and fill-in you mention.

Data generation

Since no data was provided in the question let us generate such data.

Dimensions[smat]
(* {15000, 15000} *)

smat["Density"]*Apply[Times, Dimensions[smat]]
(* 800000. *)

exportRows =
Map[{#[[1, 1, 1]],
Transpose[{#[[All, 1, 2]], RandomReal[{0, 10}, Length[#]]}]} &,
GatherBy[Most[ArrayRules[smat]], #[[1, 1]] &]];

aijText =
StringJoin@
Riffle[Map[
ToString[#[]] <> ":" <>
StringReplace[
StringTake[ToString[#[]], {2, -2}], {"{" -> "(",
"}" -> ")", "," ->" " }] &, exportRows], "\n"];

StringLength[aijText]

(* 13872182 *)

Conversion

Here is the conversion code. Note the use of Cases to eliminate wrong parsings. (I was too lazy to identify the reason, for smaller matrices the process happens without problems...)

AbsoluteTiming[
t1 = StringReplace[
aijText, {")" ~~ (WhitespaceCharacter ..) ~~ "(" -> "},{",
"(" -> "{", ")" -> "}",
x : NumberString ~~ (WhitespaceCharacter ..) ~~ y : NumberString :> x <> "," <> y}];
]

(* {1.534615, Null} *)

AbsoluteTiming[
triplets = Flatten[#, 1] &@
Map[{ToExpression[#[]],
ToExpression["{" <> #[] <> "}"]} &,
Map[StringCases[#, x__ ~~ ":" ~~ y__ :> {x, y}][] &,
StringSplit[t1, "\n"]]]];
]

(* {2.763, Null} *)

triplets // Length
(* 799893 *)

AbsoluteTiming[
triplets = Cases[triplets, {_?NumberQ, _?NumberQ, _?NumberQ}];
]

(* {0.520321, Null} *)

triplets // Length
(* 799893 *)

smat2 =
SparseArray[Thread[#[[All, 1 ;; 2]] -> #[[All, 3]]] &@triplets]
• The example in the question doesn't imply , between adjacent ) and (. Hence I would replace StringReplace[aijText, {"(" -> "{", ")" -> "}"}] with StringReplace[aijText, {") (" -> "},{", "(" -> "{", ")" -> "}"}] (assuming one space between adjacent ) and (). May 29 '17 at 14:17
• @AlexeyPopkov Thank you, for looking into this, I made the related corrections. May 29 '17 at 19:37
• @AntonAntonov I just had to add triplets = Map[{#[] + 1, #[] + 1, #[]} &, triplets] to make this work, as the original file starts with zero, and Mathematica complains with SparseArray::posr: The left-hand side of {0,0}->1 in {{0,0}->1,{0,2}->2,{1,1}->1,{2,0}->-(1/1000),{2,2}->2} is not a position or a pattern that will match the position of an element in an array with depth 2.. Thanks! May 30 '17 at 10:15
• @senseiwa Ok, I am glad my answer was useful. May 30 '17 at 13:17