<html><head><meta http-equiv="content-type" content="text/html; charset=us-ascii"></head><body style="overflow-wrap: break-word; -webkit-nbsp-mode: space; line-break: after-white-space;"><div><br></div> Are the nonzero structures of all the rows related? If they are, one could devise a routine to take advantage of this relationship, but if the nonzero structures of each row are "randomly" different from all the other rows, then it is difficult to see how one can take advantage of the sparsity.<div><br></div><div><br><div><br><blockquote type="cite"><div>On Aug 29, 2023, at 12:50 PM, Thanasis Boutsikakis <thanasis.boutsikakis@corintis.com> wrote:</div><br class="Apple-interchange-newline"><div><meta http-equiv="content-type" content="text/html; charset=us-ascii"><div style="overflow-wrap: break-word; -webkit-nbsp-mode: space; line-break: after-white-space;">Hi all, I have the following code that orthogonalizes a PETSc matrix. The problem is that this implementation requires that the PETSc matrix is dense, otherwise, it fails at bv.SetFromOptions(). Hence the assert in orthogonality().<div><br></div><div>What could I do in order to be able to orthogonalize sparse matrices as well? Could I convert it efficiently? (I tried to no avail)</div><div><br></div><div>Thanks!<br><div><br></div><div><div style="color: rgb(29, 33, 38); background-color: rgb(255, 255, 255); font-family: Menlo, Monaco, "Courier New", monospace; line-height: 18px; white-space: pre;"><div><div style="line-height: 18px;"><div><span style="color: rgb(10, 66, 131);">"""Experimenting with matrix orthogonalization"""</span></div><br><div><span style="color: rgb(168, 3, 20);">import</span> <span style="color: rgb(72, 31, 148);">contextlib</span></div><div><span style="color: rgb(168, 3, 20);">import</span> <span style="color: rgb(72, 31, 148);">sys</span></div><div><span style="color: rgb(168, 3, 20);">import</span> <span style="color: rgb(72, 31, 148);">time</span></div><div><span style="color: rgb(168, 3, 20);">import</span> <span style="color: rgb(72, 31, 148);">numpy</span> <span style="color: rgb(168, 3, 20);">as</span> <span style="color: rgb(72, 31, 148);">np</span></div><div><span style="color: rgb(168, 3, 20);">from</span> <span style="color: rgb(72, 31, 148);">firedrake</span> <span style="color: rgb(168, 3, 20);">import</span> <span style="color: rgb(13, 92, 181);">COMM_WORLD</span></div><div><span style="color: rgb(168, 3, 20);">from</span> <span style="color: rgb(72, 31, 148);">firedrake</span>.<span style="color: rgb(72, 31, 148);">petsc</span> <span style="color: rgb(168, 3, 20);">import</span> <span style="color: rgb(72, 31, 148);">PETSc</span></div><br><div><span style="color: rgb(168, 3, 20);">import</span> <span style="color: rgb(72, 31, 148);">slepc4py</span></div><br><div><span style="color: rgb(72, 31, 148);">slepc4py</span>.<span style="color: rgb(72, 31, 148);">init</span>(<span style="color: rgb(72, 31, 148);">sys</span>.argv)</div><div><span style="color: rgb(168, 3, 20);">from</span> <span style="color: rgb(72, 31, 148);">slepc4py</span> <span style="color: rgb(168, 3, 20);">import</span> <span style="color: rgb(72, 31, 148);">SLEPc</span></div><br><div><span style="color: rgb(168, 3, 20);">from</span> <span style="color: rgb(72, 31, 148);">numpy</span>.<span style="color: rgb(72, 31, 148);">testing</span> <span style="color: rgb(168, 3, 20);">import</span> <span style="color: rgb(72, 31, 148);">assert_array_almost_equal</span></div><br><div><span style="color: rgb(13, 92, 181);">EPSILON_USER</span> <span style="color: rgb(168, 3, 20);">=</span> <span style="color: rgb(13, 92, 181);">1e-4</span></div><div><span style="color: rgb(13, 92, 181);">EPS</span> <span style="color: rgb(168, 3, 20);">=</span> <span style="color: rgb(72, 31, 148);">sys</span>.float_info.epsilon</div><br><br><div><span style="color: rgb(168, 3, 20);">def</span> <span style="color: rgb(72, 31, 148);">Print</span>(<span style="color: rgb(166, 61, 0);">message</span>: <span style="color: rgb(72, 31, 148);">str</span>):</div><div> <span style="color: rgb(10, 66, 131);">"""Print function that prints only on rank 0 with color</span></div><br><div><span style="color: rgb(10, 66, 131);"> Args:</span></div><div><span style="color: rgb(10, 66, 131);"> message (str): message to be printed</span></div><div><span style="color: rgb(10, 66, 131);"> """</span></div><div> <span style="color: rgb(72, 31, 148);">PETSc</span>.Sys.Print(<span style="color: rgb(166, 61, 0);">message</span>)</div><br><br><div><span style="color: rgb(168, 3, 20);">def</span> <span style="color: rgb(72, 31, 148);">create_petsc_matrix</span>(<span style="color: rgb(166, 61, 0);">input_array</span>, <span style="color: rgb(166, 61, 0);">sparse</span><span style="color: rgb(168, 3, 20);">=</span><span style="color: rgb(13, 92, 181);">True</span>):</div><div> <span style="color: rgb(10, 66, 131);">"""Create a PETSc matrix from an input_array</span></div><br><div><span style="color: rgb(10, 66, 131);"> Args:</span></div><div><span style="color: rgb(10, 66, 131);"> input_array (np array): Input array</span></div><div><span style="color: rgb(10, 66, 131);"> partition_like (PETSc mat, optional): Petsc matrix. Defaults to None.</span></div><div><span style="color: rgb(10, 66, 131);"> sparse (bool, optional): Toggle for sparese or dense. Defaults to True.</span></div><br><div><span style="color: rgb(10, 66, 131);"> Returns:</span></div><div><span style="color: rgb(10, 66, 131);"> PETSc mat: PETSc matrix</span></div><div><span style="color: rgb(10, 66, 131);"> """</span></div><div> <span style="color: rgb(164, 176, 188);"># Check if input_array is 1D and reshape if necessary</span></div><div> <span style="color: rgb(168, 3, 20);">assert</span> <span style="color: rgb(72, 31, 148);">len</span>(<span style="color: rgb(166, 61, 0);">input_array</span>.shape) <span style="color: rgb(168, 3, 20);">==</span> <span style="color: rgb(13, 92, 181);">2</span>, <span style="color: rgb(10, 66, 131);">"Input array should be 2-dimensional"</span></div><div> global_rows, global_cols <span style="color: rgb(168, 3, 20);">=</span> <span style="color: rgb(166, 61, 0);">input_array</span>.shape</div><br><div> size <span style="color: rgb(168, 3, 20);">=</span> ((<span style="color: rgb(13, 92, 181);">None</span>, global_rows), (global_cols, global_cols))</div><br><div> <span style="color: rgb(164, 176, 188);"># Create a sparse or dense matrix based on the 'sparse' argument</span></div><div> <span style="color: rgb(168, 3, 20);">if</span> <span style="color: rgb(166, 61, 0);">sparse</span>:</div><div> matrix <span style="color: rgb(168, 3, 20);">=</span> <span style="color: rgb(72, 31, 148);">PETSc</span>.Mat().createAIJ(<span style="color: rgb(166, 61, 0);">size</span><span style="color: rgb(168, 3, 20);">=</span>size, <span style="color: rgb(166, 61, 0);">comm</span><span style="color: rgb(168, 3, 20);">=</span><span style="color: rgb(13, 92, 181);">COMM_WORLD</span>)</div><div> <span style="color: rgb(168, 3, 20);">else</span>:</div><div> matrix <span style="color: rgb(168, 3, 20);">=</span> <span style="color: rgb(72, 31, 148);">PETSc</span>.Mat().createDense(<span style="color: rgb(166, 61, 0);">size</span><span style="color: rgb(168, 3, 20);">=</span>size, <span style="color: rgb(166, 61, 0);">comm</span><span style="color: rgb(168, 3, 20);">=</span><span style="color: rgb(13, 92, 181);">COMM_WORLD</span>)</div><div> matrix.setUp()</div><br><div> local_rows_start, local_rows_end <span style="color: rgb(168, 3, 20);">=</span> matrix.getOwnershipRange()</div><br><div> <span style="color: rgb(168, 3, 20);">for</span> counter, i <span style="color: rgb(168, 3, 20);">in</span> <span style="color: rgb(72, 31, 148);">enumerate</span>(<span style="color: rgb(72, 31, 148);">range</span>(local_rows_start, local_rows_end)):</div><div> <span style="color: rgb(164, 176, 188);"># Calculate the correct row in the array for the current process</span></div><div> row_in_array <span style="color: rgb(168, 3, 20);">=</span> counter <span style="color: rgb(168, 3, 20);">+</span> local_rows_start</div><div> matrix.setValues(</div><div> i, <span style="color: rgb(72, 31, 148);">range</span>(global_cols), <span style="color: rgb(166, 61, 0);">input_array</span>[row_in_array, :], <span style="color: rgb(166, 61, 0);">addv</span><span style="color: rgb(168, 3, 20);">=</span><span style="color: rgb(13, 92, 181);">False</span></div><div> )</div><br><div> <span style="color: rgb(164, 176, 188);"># Assembly the matrix to compute the final structure</span></div><div> matrix.assemblyBegin()</div><div> matrix.assemblyEnd()</div><br><div> <span style="color: rgb(168, 3, 20);">return</span> matrix</div><br><br><div><span style="color: rgb(168, 3, 20);">def</span> <span style="color: rgb(72, 31, 148);">orthogonality</span>(<span style="color: rgb(166, 61, 0);">A</span>): <span style="color: rgb(164, 176, 188);"># sourcery skip: avoid-builtin-shadow</span></div><div> <span style="color: rgb(10, 66, 131);">"""Checking and correcting orthogonality</span></div><br><div><span style="color: rgb(10, 66, 131);"> Args:</span></div><div><span style="color: rgb(10, 66, 131);"> A (PETSc.Mat): Matrix of size [m x k].</span></div><br><div><span style="color: rgb(10, 66, 131);"> Returns:</span></div><div><span style="color: rgb(10, 66, 131);"> PETSc.Mat: Matrix of size [m x k].</span></div><div><span style="color: rgb(10, 66, 131);"> """</span></div><div> <span style="color: rgb(164, 176, 188);"># Check if the matrix is dense</span></div><div> mat_type <span style="color: rgb(168, 3, 20);">=</span> <span style="color: rgb(166, 61, 0);">A</span>.getType()</div><div> <span style="color: rgb(168, 3, 20);">assert</span> mat_type <span style="color: rgb(168, 3, 20);">in</span> (</div><div> <span style="color: rgb(10, 66, 131);">"seqdense"</span>,</div><div> <span style="color: rgb(10, 66, 131);">"mpidense"</span>,</div><div> ), <span style="color: rgb(10, 66, 131);">"A must be a dense matrix. SLEPc.BV().createFromMat() requires a dense matrix."</span></div><br><div> m, k <span style="color: rgb(168, 3, 20);">=</span> <span style="color: rgb(166, 61, 0);">A</span>.getSize()</div><br><div> Phi1 <span style="color: rgb(168, 3, 20);">=</span> <span style="color: rgb(166, 61, 0);">A</span>.getColumnVector(<span style="color: rgb(13, 92, 181);">0</span>)</div><div> Phi2 <span style="color: rgb(168, 3, 20);">=</span> <span style="color: rgb(166, 61, 0);">A</span>.getColumnVector(k <span style="color: rgb(168, 3, 20);">-</span> <span style="color: rgb(13, 92, 181);">1</span>)</div><br><div> <span style="color: rgb(164, 176, 188);"># Compute dot product using PETSc function</span></div><div> dot_product <span style="color: rgb(168, 3, 20);">=</span> Phi1.dot(Phi2)</div><br><div> <span style="color: rgb(168, 3, 20);">if</span> <span style="color: rgb(72, 31, 148);">abs</span>(dot_product) <span style="color: rgb(168, 3, 20);">></span> <span style="color: rgb(72, 31, 148);">min</span>(<span style="color: rgb(13, 92, 181);">EPSILON_USER</span>, <span style="color: rgb(13, 92, 181);">EPS</span> <span style="color: rgb(168, 3, 20);">*</span> m):</div><div> <span style="color: rgb(72, 31, 148);">Print</span>(<span style="color: rgb(10, 66, 131);">" Matrix is not orthogonal"</span>)</div><br><div> <span style="color: rgb(164, 176, 188);"># Type can be CHOL, GS, mro(), SVQB, TSQR, TSQRCHOL</span></div><div> _type <span style="color: rgb(168, 3, 20);">=</span> <span style="color: rgb(72, 31, 148);">SLEPc</span>.BV().OrthogBlockType.<span style="color: rgb(13, 92, 181);">GS</span></div><br><div> bv <span style="color: rgb(168, 3, 20);">=</span> <span style="color: rgb(72, 31, 148);">SLEPc</span>.BV().createFromMat(<span style="color: rgb(166, 61, 0);">A</span>)</div><div> bv.setFromOptions()</div><div> bv.setOrthogonalization(_type)</div><div> bv.orthogonalize()</div><br><div> <span style="color: rgb(166, 61, 0);">A</span> <span style="color: rgb(168, 3, 20);">=</span> bv.createMat()</div><br><div> <span style="color: rgb(72, 31, 148);">Print</span>(<span style="color: rgb(10, 66, 131);">" Matrix successfully orthogonalized"</span>)</div><br><div> <span style="color: rgb(164, 176, 188);"># # Assembly the matrix to compute the final structure</span></div><div> <span style="color: rgb(168, 3, 20);">if</span> <span style="color: rgb(168, 3, 20);">not</span> <span style="color: rgb(166, 61, 0);">A</span>.assembled:</div><div> <span style="color: rgb(166, 61, 0);">A</span>.assemblyBegin()</div><div> <span style="color: rgb(166, 61, 0);">A</span>.assemblyEnd()</div><div> <span style="color: rgb(168, 3, 20);">else</span>:</div><div> <span style="color: rgb(72, 31, 148);">Print</span>(<span style="color: rgb(10, 66, 131);">" Matrix is orthogonal"</span>)</div><br><div> <span style="color: rgb(168, 3, 20);">return</span> <span style="color: rgb(166, 61, 0);">A</span></div><br><br><div><span style="color: rgb(164, 176, 188);"># --------------------------------------------</span></div><div><span style="color: rgb(164, 176, 188);"># EXP: Orthogonalization of an mpi PETSc matrix</span></div><div><span style="color: rgb(164, 176, 188);"># --------------------------------------------</span></div><br><div>m, k <span style="color: rgb(168, 3, 20);">=</span> <span style="color: rgb(13, 92, 181);">11</span>, <span style="color: rgb(13, 92, 181);">7</span></div><div><span style="color: rgb(164, 176, 188);"># Generate the random numpy matrices</span></div><div><span style="color: rgb(72, 31, 148);">np</span>.<span style="color: rgb(72, 31, 148);">random</span>.seed(<span style="color: rgb(13, 92, 181);">0</span>) <span style="color: rgb(164, 176, 188);"># sets the seed to 0</span></div><div>A_np <span style="color: rgb(168, 3, 20);">=</span> <span style="color: rgb(72, 31, 148);">np</span>.<span style="color: rgb(72, 31, 148);">random</span>.randint(<span style="color: rgb(166, 61, 0);">low</span><span style="color: rgb(168, 3, 20);">=</span><span style="color: rgb(13, 92, 181);">0</span>, <span style="color: rgb(166, 61, 0);">high</span><span style="color: rgb(168, 3, 20);">=</span><span style="color: rgb(13, 92, 181);">6</span>, <span style="color: rgb(166, 61, 0);">size</span><span style="color: rgb(168, 3, 20);">=</span>(m, k))</div><br><div><span style="color: rgb(13, 92, 181);">A</span> <span style="color: rgb(168, 3, 20);">=</span> <span style="color: rgb(72, 31, 148);">create_petsc_matrix</span>(A_np, <span style="color: rgb(166, 61, 0);">sparse</span><span style="color: rgb(168, 3, 20);">=</span><span style="color: rgb(13, 92, 181);">False</span>)</div><br><div>A_orthogonal <span style="color: rgb(168, 3, 20);">=</span> <span style="color: rgb(72, 31, 148);">orthogonality</span>(<span style="color: rgb(13, 92, 181);">A</span>)</div><br><div><span style="color: rgb(164, 176, 188);"># --------------------------------------------</span></div><div><span style="color: rgb(164, 176, 188);"># TEST: Orthogonalization of a numpy matrix</span></div><div><span style="color: rgb(164, 176, 188);"># --------------------------------------------</span></div><div><span style="color: rgb(164, 176, 188);"># Generate A_np_orthogonal</span></div><div>A_np_orthogonal, _ <span style="color: rgb(168, 3, 20);">=</span> <span style="color: rgb(72, 31, 148);">np</span>.<span style="color: rgb(72, 31, 148);">linalg</span>.<span style="color: rgb(72, 31, 148);">qr</span>(A_np)</div><br><div><span style="color: rgb(164, 176, 188);"># Get the local values from A_orthogonal</span></div><div>local_rows_start, local_rows_end <span style="color: rgb(168, 3, 20);">=</span> A_orthogonal.getOwnershipRange()</div><div>A_orthogonal_local <span style="color: rgb(168, 3, 20);">=</span> A_orthogonal.getValues(</div><div> <span style="color: rgb(72, 31, 148);">range</span>(local_rows_start, local_rows_end), <span style="color: rgb(72, 31, 148);">range</span>(k)</div><div>)</div><br><div><span style="color: rgb(164, 176, 188);"># Assert the correctness of the multiplication for the local subset</span></div><div><span style="color: rgb(72, 31, 148);">assert_array_almost_equal</span>(</div><div> <span style="color: rgb(72, 31, 148);">np</span>.abs(A_orthogonal_local),</div><div> <span style="color: rgb(72, 31, 148);">np</span>.abs(A_np_orthogonal[local_rows_start:local_rows_end, :]),</div><div> <span style="color: rgb(166, 61, 0);">decimal</span><span style="color: rgb(168, 3, 20);">=</span><span style="color: rgb(13, 92, 181);">5</span>,</div><div>)</div></div></div></div></div></div></div></div></blockquote></div><br></div></body></html>