ceed-basis.c (2efebffe38fa9227caaeab504b43e3a698cb86d7) ceed-basis.c (2c2ea1dbee80fceecd2c97f30b09f8c87820a53e)
1// Copyright (c) 2017-2022, Lawrence Livermore National Security, LLC and other CEED contributors.
2// All Rights Reserved. See the top-level LICENSE and NOTICE files for details.
3//
4// SPDX-License-Identifier: BSD-2-Clause
5//
6// This file is part of CEED: http://github.com/ceed
7
8#include <ceed-impl.h>

--- 660 unchanged lines hidden (view full) ---

669 @param[in,out] mat Row-major matrix to be factorized in place
670 @param[out] lambda Vector of length n of eigenvalues
671 @param[in] n Number of rows/columns
672
673 @return An error code: 0 - success, otherwise - failure
674
675 @ref Utility
676**/
1// Copyright (c) 2017-2022, Lawrence Livermore National Security, LLC and other CEED contributors.
2// All Rights Reserved. See the top-level LICENSE and NOTICE files for details.
3//
4// SPDX-License-Identifier: BSD-2-Clause
5//
6// This file is part of CEED: http://github.com/ceed
7
8#include <ceed-impl.h>

--- 660 unchanged lines hidden (view full) ---

669 @param[in,out] mat Row-major matrix to be factorized in place
670 @param[out] lambda Vector of length n of eigenvalues
671 @param[in] n Number of rows/columns
672
673 @return An error code: 0 - success, otherwise - failure
674
675 @ref Utility
676**/
677CeedPragmaOptimizeOff int CeedSymmetricSchurDecomposition(Ceed ceed, CeedScalar *mat, CeedScalar *lambda, CeedInt n) {
677CeedPragmaOptimizeOff
678int CeedSymmetricSchurDecomposition(Ceed ceed, CeedScalar *mat, CeedScalar *lambda, CeedInt n) {
678 // Check bounds for clang-tidy
679 if (n < 2) {
680 // LCOV_EXCL_START
681 return CeedError(ceed, CEED_ERROR_UNSUPPORTED, "Cannot compute symmetric Schur decomposition of scalars");
682 // LCOV_EXCL_STOP
683 }
684
685 CeedScalar v[n - 1], tau[n - 1], mat_T[n * n];

--- 109 unchanged lines hidden (view full) ---

795 // Check convergence
796 if (itr == max_itr && q < n - 1) {
797 // LCOV_EXCL_START
798 return CeedError(ceed, CEED_ERROR_MINOR, "Symmetric QR failed to converge");
799 // LCOV_EXCL_STOP
800 }
801 return CEED_ERROR_SUCCESS;
802}
679 // Check bounds for clang-tidy
680 if (n < 2) {
681 // LCOV_EXCL_START
682 return CeedError(ceed, CEED_ERROR_UNSUPPORTED, "Cannot compute symmetric Schur decomposition of scalars");
683 // LCOV_EXCL_STOP
684 }
685
686 CeedScalar v[n - 1], tau[n - 1], mat_T[n * n];

--- 109 unchanged lines hidden (view full) ---

796 // Check convergence
797 if (itr == max_itr && q < n - 1) {
798 // LCOV_EXCL_START
799 return CeedError(ceed, CEED_ERROR_MINOR, "Symmetric QR failed to converge");
800 // LCOV_EXCL_STOP
801 }
802 return CEED_ERROR_SUCCESS;
803}
803CeedPragmaOptimizeOn;
804CeedPragmaOptimizeOn
804
805/**
806 @brief Return Simultaneous Diagonalization of two matrices.
807
808 This solves the generalized eigenvalue problem A x = lambda B x, where A and B are symmetric and B is positive definite.
809 We generate the matrix X and vector Lambda such that X^T A X = Lambda and X^T B X = I.
810 This is equivalent to the LAPACK routine 'sygv' with TYPE = 1.
811
812 @param[in] ceed Ceed context for error handling
813 @param[in] mat_A Row-major matrix to be factorized with eigenvalues
814 @param[in] mat_B Row-major matrix to be factorized to identity
815 @param[out] mat_X Row-major orthogonal matrix
816 @param[out] lambda Vector of length n of generalized eigenvalues
817 @param[in] n Number of rows/columns
818
819 @return An error code: 0 - success, otherwise - failure
820
821 @ref Utility
822**/
805
806/**
807 @brief Return Simultaneous Diagonalization of two matrices.
808
809 This solves the generalized eigenvalue problem A x = lambda B x, where A and B are symmetric and B is positive definite.
810 We generate the matrix X and vector Lambda such that X^T A X = Lambda and X^T B X = I.
811 This is equivalent to the LAPACK routine 'sygv' with TYPE = 1.
812
813 @param[in] ceed Ceed context for error handling
814 @param[in] mat_A Row-major matrix to be factorized with eigenvalues
815 @param[in] mat_B Row-major matrix to be factorized to identity
816 @param[out] mat_X Row-major orthogonal matrix
817 @param[out] lambda Vector of length n of generalized eigenvalues
818 @param[in] n Number of rows/columns
819
820 @return An error code: 0 - success, otherwise - failure
821
822 @ref Utility
823**/
823CeedPragmaOptimizeOff int CeedSimultaneousDiagonalization(Ceed ceed, CeedScalar *mat_A, CeedScalar *mat_B, CeedScalar *mat_X, CeedScalar *lambda,
824 CeedInt n) {
824CeedPragmaOptimizeOff
825int CeedSimultaneousDiagonalization(Ceed ceed, CeedScalar *mat_A, CeedScalar *mat_B, CeedScalar *mat_X, CeedScalar *lambda, CeedInt n) {
825 CeedScalar *mat_C, *mat_G, *vec_D;
826 CeedCall(CeedCalloc(n * n, &mat_C));
827 CeedCall(CeedCalloc(n * n, &mat_G));
828 CeedCall(CeedCalloc(n, &vec_D));
829
830 // Compute B = G D G^T
831 memcpy(mat_G, mat_B, n * n * sizeof(mat_B[0]));
832 CeedCall(CeedSymmetricSchurDecomposition(ceed, mat_G, vec_D, n));

--- 57 unchanged lines hidden (view full) ---

890 CeedCall(CeedMatrixMatrixMultiply(ceed, (const CeedScalar *)mat_G, (const CeedScalar *)mat_C, mat_X, n, n, n));
891
892 // Cleanup
893 CeedCall(CeedFree(&mat_C));
894 CeedCall(CeedFree(&mat_G));
895 CeedCall(CeedFree(&vec_D));
896 return CEED_ERROR_SUCCESS;
897}
826 CeedScalar *mat_C, *mat_G, *vec_D;
827 CeedCall(CeedCalloc(n * n, &mat_C));
828 CeedCall(CeedCalloc(n * n, &mat_G));
829 CeedCall(CeedCalloc(n, &vec_D));
830
831 // Compute B = G D G^T
832 memcpy(mat_G, mat_B, n * n * sizeof(mat_B[0]));
833 CeedCall(CeedSymmetricSchurDecomposition(ceed, mat_G, vec_D, n));

--- 57 unchanged lines hidden (view full) ---

891 CeedCall(CeedMatrixMatrixMultiply(ceed, (const CeedScalar *)mat_G, (const CeedScalar *)mat_C, mat_X, n, n, n));
892
893 // Cleanup
894 CeedCall(CeedFree(&mat_C));
895 CeedCall(CeedFree(&mat_G));
896 CeedCall(CeedFree(&vec_D));
897 return CEED_ERROR_SUCCESS;
898}
898CeedPragmaOptimizeOn;
899CeedPragmaOptimizeOn
899
900/// @}
901
902/// ----------------------------------------------------------------------------
903/// CeedBasis Public API
904/// ----------------------------------------------------------------------------
905/// @addtogroup CeedBasisUser
906/// @{

--- 1110 unchanged lines hidden --- 		900
901/// @}
902
903/// ----------------------------------------------------------------------------
904/// CeedBasis Public API
905/// ----------------------------------------------------------------------------
906/// @addtogroup CeedBasisUser
907/// @{

--- 1110 unchanged lines hidden ---