14411cf47Sjeremylt /// @file
252bfb9bbSJeremy L Thompson /// Test QR Factorization
352bfb9bbSJeremy L Thompson /// \test Test QR Factorization
4*f85e4a7bSJeremy L Thompson
5*f85e4a7bSJeremy L Thompson //TESTARGS(only="cpu") {ceed_resource}
657c64913Sjeremylt #include <ceed.h>
79ac7b42eSJeremy L Thompson #include <ceed/backend.h>
89ac7b42eSJeremy L Thompson #include <math.h>
949aac155SJeremy L Thompson #include <stdio.h>
1057c64913Sjeremylt
main(int argc,char ** argv)1157c64913Sjeremylt int main(int argc, char **argv) {
1257c64913Sjeremylt Ceed ceed;
139ac7b42eSJeremy L Thompson CeedScalar A[12] = {1, -1, 4, 1, 4, -2, 1, 4, 2, 1, -1, 0};
1452bfb9bbSJeremy L Thompson CeedScalar qr[12] = {1, -1, 4, 1, 4, -2, 1, 4, 2, 1, -1, 0};
159ac7b42eSJeremy L Thompson CeedScalar A_qr[12] = {0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0};
169ac7b42eSJeremy L Thompson CeedScalar tau[4];
1757c64913Sjeremylt
1857c64913Sjeremylt CeedInit(argv[1], &ceed);
19aedaa0e5Sjeremylt
2052bfb9bbSJeremy L Thompson CeedQRFactorization(ceed, qr, tau, 4, 3);
212b730f8bSJeremy L Thompson for (CeedInt i = 0; i < 3; i++) {
222b730f8bSJeremy L Thompson for (CeedInt j = i; j < 3; j++) A_qr[i * 3 + j] = qr[i * 3 + j];
232b730f8bSJeremy L Thompson }
249ac7b42eSJeremy L Thompson CeedHouseholderApplyQ(A_qr, qr, tau, CEED_NOTRANSPOSE, 4, 3, 3, 3, 1);
259ac7b42eSJeremy L Thompson
262b730f8bSJeremy L Thompson for (CeedInt i = 0; i < 12; i++) {
272b730f8bSJeremy L Thompson if (fabs(A_qr[i] - A[i]) > 100. * CEED_EPSILON) {
289ac7b42eSJeremy L Thompson // LCOV_EXCL_START
292b730f8bSJeremy L Thompson printf("Error in QR factorization A_qr[%" CeedInt_FMT "] = %f != A[%" CeedInt_FMT "] = %f\n", i, A_qr[i], i, A[i]);
309ac7b42eSJeremy L Thompson // LCOV_EXCL_STOP
312b730f8bSJeremy L Thompson }
322b730f8bSJeremy L Thompson }
339ac7b42eSJeremy L Thompson
3457c64913Sjeremylt CeedDestroy(&ceed);
3557c64913Sjeremylt return 0;
3657c64913Sjeremylt }
37