used CeedMatrixPseudoinverse in CeedBasisGetCollocatedGrad didn't work

interface/ceed-basis.c

@@ -258,37 +258,25 @@ static int CeedBasisCreateProjectionMatrices(CeedBasis basis_from, CeedBasis bas
   @ref Backend
 **/
 int CeedBasisGetCollocatedGrad(CeedBasis basis, CeedScalar *collo_grad_1d) {
-  int         i, j, k;
   Ceed        ceed;
   CeedInt     P_1d = (basis)->P_1d, Q_1d = (basis)->Q_1d;
-  CeedScalar *interp_1d, *grad_1d, *tau;
+  CeedScalar *interp_1d, *interp_1d_pinv, *grad_1d, *tau;
 
   CeedCall(CeedMalloc(Q_1d * P_1d, &interp_1d));
+  CeedCall(CeedMalloc(P_1d * Q_1d, &interp_1d_pinv));
   CeedCall(CeedMalloc(Q_1d * P_1d, &grad_1d));
   CeedCall(CeedMalloc(Q_1d, &tau));
   memcpy(interp_1d, (basis)->interp_1d, Q_1d * P_1d * sizeof(basis)->interp_1d[0]);
   memcpy(grad_1d, (basis)->grad_1d, Q_1d * P_1d * sizeof(basis)->interp_1d[0]);
 
   // QR Factorization, interp_1d = Q R
   CeedCall(CeedBasisGetCeed(basis, &ceed));
-  CeedCall(CeedQRFactorization(ceed, interp_1d, tau, Q_1d, P_1d));
-  // Note: This function is for backend use, so all errors are terminal and we do not need to clean up memory on failure.
-
-  // Apply Rinv, collo_grad_1d = grad_1d Rinv
-  for (i = 0; i < Q_1d; i++) {  // Row i
-    collo_grad_1d[Q_1d * i] = grad_1d[P_1d * i] / interp_1d[0];
-    for (j = 1; j < P_1d; j++) {  // Column j
-      collo_grad_1d[j + Q_1d * i] = grad_1d[j + P_1d * i];
-      for (k = 0; k < j; k++) collo_grad_1d[j + Q_1d * i] -= interp_1d[j + P_1d * k] * collo_grad_1d[k + Q_1d * i];
-      collo_grad_1d[j + Q_1d * i] /= interp_1d[j + P_1d * j];
-    }
-    for (j = P_1d; j < Q_1d; j++) collo_grad_1d[j + Q_1d * i] = 0;
-  }
 
-  // Apply Q^T, collo_grad_1d = collo_grad_1d Q^T
-  CeedCall(CeedHouseholderApplyQ(collo_grad_1d, interp_1d, tau, CEED_NOTRANSPOSE, Q_1d, Q_1d, P_1d, 1, Q_1d));
+  CeedCall(CeedMatrixPseudoinverse(ceed, interp_1d, Q_1d, P_1d, interp_1d_pinv));
+  CeedCall(CeedMatrixMatrixMultiply(ceed, (const CeedScalar *)interp_1d_pinv, (const CeedScalar *)grad_1d, collo_grad_1d, Q_1d, P_1d, P_1d));
 
   CeedCall(CeedFree(&interp_1d));
+  CeedCall(CeedFree(&interp_1d_pinv));
   CeedCall(CeedFree(&grad_1d));
   CeedCall(CeedFree(&tau));
   return CEED_ERROR_SUCCESS;
@@ -652,24 +640,23 @@ int CeedHouseholderApplyQ(CeedScalar *mat_A, const CeedScalar *mat_Q, const Ceed
 int CeedMatrixPseudoinverse(Ceed ceed, CeedScalar *mat, CeedInt m, CeedInt n, CeedScalar *mat_pinv) {
   CeedScalar *tau, *I;
 
-  CeedCall(CeedCalloc(m * n, &I));
+  CeedCall(CeedCalloc(m * m, &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;
+  for (CeedInt i = 0; i < m; i++) I[i * m + i] = 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];
+    mat_pinv[m * i] = I[n * 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] = I[j + n * i];
+      for (CeedInt k = 0; k < j; k++) mat_pinv[j + m * i] -= mat[j + n * 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));
+  CeedCall(CeedHouseholderApplyQ(mat, mat_pinv, tau, CEED_NOTRANSPOSE, m, n, n, 1, n));
 
   return CEED_ERROR_SUCCESS;
 }