Added CeedMatrixPseudoinverse utility

CEED · Jul 10, 2023 · cc3a690 · cc3a690
1 parent 63ca180
commit cc3a690
Show file tree

Hide file tree

Showing 2 changed files with 39 additions and 0 deletions.
diff --git a/include/ceed/backend.h b/include/ceed/backend.h
@@ -341,6 +341,7 @@ CEED_INTERN int CeedMatrixMatrixMultiply(Ceed ceed, const CeedScalar *mat_A, con
 CEED_EXTERN int CeedQRFactorization(Ceed ceed, CeedScalar *mat, CeedScalar *tau, CeedInt m, CeedInt n);
 CEED_EXTERN int CeedHouseholderApplyQ(CeedScalar *mat_A, const CeedScalar *mat_Q, const CeedScalar *tau, CeedTransposeMode t_mode, CeedInt m,
                                       CeedInt n, CeedInt k, CeedInt row, CeedInt col);
+CEED_EXTERN int CeedMatrixPseudoinverse(Ceed ceed, CeedScalar *mat, CeedInt m, CeedInt n, CeedScalar *mat_pinv);
 CEED_EXTERN int CeedSymmetricSchurDecomposition(Ceed ceed, CeedScalar *mat, CeedScalar *lambda, CeedInt n);
 CEED_EXTERN int CeedSimultaneousDiagonalization(Ceed ceed, CeedScalar *mat_A, CeedScalar *mat_B, CeedScalar *x, CeedScalar *lambda, CeedInt n);
 

diff --git a/interface/ceed-basis.c b/interface/ceed-basis.c
@@ -636,6 +636,44 @@ int CeedHouseholderApplyQ(CeedScalar *mat_A, const CeedScalar *mat_Q, const Ceed
   return CEED_ERROR_SUCCESS;
 }
 
+/**
+  @brief Return pseudoinverse of a matrix
+
+  @param[in]     ceed Ceed context for error handling
+  @param[in]     mat  Row-major matrix to be factorized in place
+  @param[in]     m    Number of rows
+  @param[in]     n    Number of columns
+  @param[out]    mat_pinv  Row-major pseudoinverse matrix
+
+  @return An error code: 0 - success, otherwise - failure
+
+  @ref Utility
+**/
+int CeedMatrixPseudoinverse(Ceed ceed, CeedScalar *mat, CeedInt m, CeedInt n, CeedScalar *mat_pinv) {
+  CeedScalar *tau, *I;
+
+  CeedCall(CeedCalloc(m * n, &I));
+  CeedCall(CeedCalloc(m, &tau));
+  // -- QR Factorization, mat = Q R
+  CeedCall(CeedQRFactorization(ceed, mat, tau, m, n));
+  // -- mat_pinv = R_inv Q^T
+  for (CeedInt i = 0; i < m; i++)
+    for (CeedInt j = 0; j < n; j++) I[i * m + j] = 1.0;
+  // ---- Apply R_inv, mat_pinv = I R_inv
+  for (CeedInt i = 0; i < m; i++) {  // Row i
+    mat_pinv[m * i] = I[m * i] / mat[0];
+    for (CeedInt j = 1; j < n; j++) {  // Column j
+      mat_pinv[j + m * i] = I[j + m * i];
+      for (CeedInt k = 0; k < j; k++) mat_pinv[j + m * i] -= mat[j + m * k] * mat_pinv[k + m * i];
+      mat_pinv[j + m * i] /= mat[j + n * j];
+    }
+  }
+  // ---- Apply Q^T, mat_pinv = R_inv Q^T
+  CeedCall(CeedHouseholderApplyQ(mat_pinv, mat, tau, CEED_NOTRANSPOSE, m, n, n, 1, m));
+
+  return CEED_ERROR_SUCCESS;
+}
+
 /**
   @brief Return symmetric Schur decomposition of the symmetric matrix mat via symmetric QR factorization