@@ -99,7 +99,8 @@ int CeedHouseholderApplyQ(CeedScalar *A, const CeedScalar *Q,
9999 CeedInt m , CeedInt n , CeedInt k ,
100100 CeedInt row , CeedInt col ) {
101101 int ierr ;
102- CeedScalar v [m ];
102+ CeedScalar * v ;
103+ ierr = CeedMalloc (m , & v ); CeedChk (ierr );
103104 for (CeedInt ii = 0 ; ii < k ; ii ++ ) {
104105 CeedInt i = t_mode == CEED_TRANSPOSE ? ii : k - 1 - ii ;
105106 for (CeedInt j = i + 1 ; j < m ; j ++ )
@@ -108,6 +109,7 @@ int CeedHouseholderApplyQ(CeedScalar *A, const CeedScalar *Q,
108109 ierr = CeedHouseholderReflect (& A [i * row ], & v [i ], tau [i ], m - i , n , row , col );
109110 CeedChk (ierr );
110111 }
112+ ierr = CeedFree (& v ); CeedChk (ierr );
111113 return CEED_ERROR_SUCCESS ;
112114}
113115
@@ -200,10 +202,11 @@ int CeedBasisGetCollocatedGrad(CeedBasis basis, CeedScalar *collo_grad_1d) {
200202 int i , j , k ;
201203 Ceed ceed ;
202204 CeedInt ierr , P_1d = (basis )-> P_1d , Q_1d = (basis )-> Q_1d ;
203- CeedScalar * interp_1d , * grad_1d , tau [ Q_1d ] ;
205+ CeedScalar * interp_1d , * grad_1d , * tau ;
204206
205207 ierr = CeedMalloc (Q_1d * P_1d , & interp_1d ); CeedChk (ierr );
206208 ierr = CeedMalloc (Q_1d * P_1d , & grad_1d ); CeedChk (ierr );
209+ ierr = CeedMalloc (Q_1d , & tau ); CeedChk (ierr );
207210 memcpy (interp_1d , (basis )-> interp_1d , Q_1d * P_1d * sizeof (basis )-> interp_1d [0 ]);
208211 memcpy (grad_1d , (basis )-> grad_1d , Q_1d * P_1d * sizeof (basis )-> interp_1d [0 ]);
209212
@@ -232,6 +235,7 @@ int CeedBasisGetCollocatedGrad(CeedBasis basis, CeedScalar *collo_grad_1d) {
232235
233236 ierr = CeedFree (& interp_1d ); CeedChk (ierr );
234237 ierr = CeedFree (& grad_1d ); CeedChk (ierr );
238+ ierr = CeedFree (& tau ); CeedChk (ierr );
235239 return CEED_ERROR_SUCCESS ;
236240}
237241
@@ -1435,37 +1439,45 @@ int CeedSimultaneousDiagonalization(Ceed ceed, CeedScalar *mat_A,
14351439 CeedScalar * mat_B , CeedScalar * x ,
14361440 CeedScalar * lambda , CeedInt n ) {
14371441 int ierr ;
1438- CeedScalar mat_C [n * n ], matG [n * n ], vecD [n ];
1442+ CeedScalar * mat_C , * mat_G , * vecD ;
1443+ ierr = CeedCalloc (n * n , & mat_C ); CeedChk (ierr );
1444+ ierr = CeedCalloc (n * n , & mat_G ); CeedChk (ierr );
1445+ ierr = CeedCalloc (n , & vecD ); CeedChk (ierr );
14391446
14401447 // Compute B = G D G^T
1441- memcpy (matG , mat_B , n * n * sizeof (mat_B [0 ]));
1442- ierr = CeedSymmetricSchurDecomposition (ceed , matG , vecD , n ); CeedChk (ierr );
1448+ memcpy (mat_G , mat_B , n * n * sizeof (mat_B [0 ]));
1449+ ierr = CeedSymmetricSchurDecomposition (ceed , mat_G , vecD , n ); CeedChk (ierr );
14431450 for (CeedInt i = 0 ; i < n ; i ++ )
14441451 vecD [i ] = sqrt (vecD [i ]);
14451452
14461453 // Compute C = (G D^1/2)^-1 A (G D^1/2)^-T
14471454 // = D^-1/2 G^T A G D^-1/2
14481455 for (CeedInt i = 0 ; i < n ; i ++ )
14491456 for (CeedInt j = 0 ; j < n ; j ++ )
1450- mat_C [j + i * n ] = matG [i + j * n ] / vecD [i ];
1457+ mat_C [j + i * n ] = mat_G [i + j * n ] / vecD [i ];
14511458 ierr = CeedMatrixMultiply (ceed , (const CeedScalar * )mat_C ,
14521459 (const CeedScalar * )mat_A , x , n , n , n );
14531460 CeedChk (ierr );
14541461 for (CeedInt i = 0 ; i < n ; i ++ )
14551462 for (CeedInt j = 0 ; j < n ; j ++ )
1456- matG [j + i * n ] = matG [j + i * n ] / vecD [j ];
1463+ mat_G [j + i * n ] = mat_G [j + i * n ] / vecD [j ];
14571464 ierr = CeedMatrixMultiply (ceed , (const CeedScalar * )x ,
1458- (const CeedScalar * )matG , mat_C , n , n , n );
1465+ (const CeedScalar * )mat_G , mat_C , n , n , n );
14591466 CeedChk (ierr );
14601467
14611468 // Compute Q^T C Q = lambda
14621469 ierr = CeedSymmetricSchurDecomposition (ceed , mat_C , lambda , n ); CeedChk (ierr );
14631470
14641471 // Set x = (G D^1/2)^-T Q
14651472 // = G D^-1/2 Q
1466- ierr = CeedMatrixMultiply (ceed , (const CeedScalar * )matG ,
1473+ ierr = CeedMatrixMultiply (ceed , (const CeedScalar * )mat_G ,
14671474 (const CeedScalar * )mat_C , x , n , n , n );
14681475 CeedChk (ierr );
1476+
1477+ // Cleanup
1478+ ierr = CeedFree (& mat_C ); CeedChk (ierr );
1479+ ierr = CeedFree (& mat_G ); CeedChk (ierr );
1480+ ierr = CeedFree (& vecD ); CeedChk (ierr );
14691481 return CEED_ERROR_SUCCESS ;
14701482}
14711483
0 commit comments