xref: /libCEED/tests/t301-basis.c (revision 49aac155e7a09736f56fb3abac0f57dab29f7cbf)
14411cf47Sjeremylt /// @file
252bfb9bbSJeremy L Thompson /// Test QR Factorization
352bfb9bbSJeremy L Thompson /// \test Test QR Factorization
457c64913Sjeremylt #include <ceed.h>
59ac7b42eSJeremy L Thompson #include <ceed/backend.h>
69ac7b42eSJeremy L Thompson #include <math.h>
7*49aac155SJeremy L Thompson #include <stdio.h>
857c64913Sjeremylt 
957c64913Sjeremylt int main(int argc, char **argv) {
1057c64913Sjeremylt   Ceed       ceed;
119ac7b42eSJeremy L Thompson   CeedScalar A[12]    = {1, -1, 4, 1, 4, -2, 1, 4, 2, 1, -1, 0};
1252bfb9bbSJeremy L Thompson   CeedScalar qr[12]   = {1, -1, 4, 1, 4, -2, 1, 4, 2, 1, -1, 0};
139ac7b42eSJeremy L Thompson   CeedScalar A_qr[12] = {0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0};
149ac7b42eSJeremy L Thompson   CeedScalar tau[4];
1557c64913Sjeremylt 
1657c64913Sjeremylt   CeedInit(argv[1], &ceed);
17aedaa0e5Sjeremylt 
1852bfb9bbSJeremy L Thompson   CeedQRFactorization(ceed, qr, tau, 4, 3);
192b730f8bSJeremy L Thompson   for (CeedInt i = 0; i < 3; i++) {
202b730f8bSJeremy L Thompson     for (CeedInt j = i; j < 3; j++) A_qr[i * 3 + j] = qr[i * 3 + j];
212b730f8bSJeremy L Thompson   }
229ac7b42eSJeremy L Thompson   CeedHouseholderApplyQ(A_qr, qr, tau, CEED_NOTRANSPOSE, 4, 3, 3, 3, 1);
239ac7b42eSJeremy L Thompson 
242b730f8bSJeremy L Thompson   for (CeedInt i = 0; i < 12; i++) {
252b730f8bSJeremy L Thompson     if (fabs(A_qr[i] - A[i]) > 100. * CEED_EPSILON) {
269ac7b42eSJeremy L Thompson       // LCOV_EXCL_START
272b730f8bSJeremy L Thompson       printf("Error in QR factorization A_qr[%" CeedInt_FMT "] = %f != A[%" CeedInt_FMT "] = %f\n", i, A_qr[i], i, A[i]);
289ac7b42eSJeremy L Thompson       // LCOV_EXCL_STOP
292b730f8bSJeremy L Thompson     }
302b730f8bSJeremy L Thompson   }
319ac7b42eSJeremy L Thompson 
3257c64913Sjeremylt   CeedDestroy(&ceed);
3357c64913Sjeremylt   return 0;
3457c64913Sjeremylt }
35