<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;">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></body></html>