diff --git a/src/MathNet.Numerics.5.1.ReSharper b/src/MathNet.Numerics.5.1.ReSharper index 1323a4a6..5b2bfd06 100644 --- a/src/MathNet.Numerics.5.1.ReSharper +++ b/src/MathNet.Numerics.5.1.ReSharper @@ -83,7 +83,9 @@ kronecker blocksize ipiv dll -Silverlight +Silverlight +namespace +da diff --git a/src/NativeWrappers/MKL/lapack.cpp b/src/NativeWrappers/MKL/lapack.cpp index 4f4686dc..4bf96054 100644 --- a/src/NativeWrappers/MKL/lapack.cpp +++ b/src/NativeWrappers/MKL/lapack.cpp @@ -1,560 +1,573 @@ #include "common.h" #include "blas.h" #include "mkl_lapack.h" +#include + extern "C" { - DLLEXPORT int s_cholesky_factor(int n, float a[]){ - char uplo = 'L'; - int info = 0; - SPOTRF(&uplo, &n, a, &n, &info); - for (int i = 0; i < n; ++i) - { - int index = i * n; - for (int j = 0; j < n && i > j; ++j) - { - a[index + j] = 0; - } - } - return info; - } - - DLLEXPORT int d_cholesky_factor(int n, double* a){ - char uplo = 'L'; - int info = 0; - DPOTRF(&uplo, &n, a, &n, &info); - for (int i = 0; i < n; ++i) - { - int index = i * n; - for (int j = 0; j < n && i > j; ++j) - { - a[index + j] = 0; - } - } - return info; - } - - DLLEXPORT int c_cholesky_factor(int n, Complex8 a[]){ - char uplo = 'L'; - int info = 0; - Complex8 zero; - zero.real = 0.0; - zero.real = 0.0; - CPOTRF(&uplo, &n, a, &n, &info); - for (int i = 0; i < n; ++i) - { - int index = i * n; - for (int j = 0; j < n && i > j; ++j) - { - a[index + j] = zero; - } - } - return info; - } - - DLLEXPORT int z_cholesky_factor(int n, Complex16 a[]){ - char uplo = 'L'; - int info = 0; - Complex16 zero; - zero.real = 0.0; - zero.real = 0.0; - ZPOTRF(&uplo, &n, a, &n, &info); - for (int i = 0; i < n; ++i) - { - int index = i * n; - for (int j = 0; j < n && i > j; ++j) - { - a[index + j] = zero; - } - } - return info; - } - - DLLEXPORT float s_matrix_norm(char norm, int m, int n, float a[], float work[]) - { - return SLANGE(&norm, &m, &n, a, &m, work); - } - - DLLEXPORT double d_matrix_norm(char norm, int m, int n, double a[], double work[]) - { - return DLANGE(&norm, &m, &n, a, &m, work); - } - - DLLEXPORT float c_matrix_norm(char norm, int m, int n, MKL_Complex8 a[], float work[]) - { - return CLANGE(&norm, &m, &n, a, &m, work); - } - - DLLEXPORT double z_matrix_norm(char norm, int m, int n, MKL_Complex16 a[], double work[]) - { - return ZLANGE(&norm, &m, &n, a, &m, work); - } - - DLLEXPORT void s_lu_factor(int m, float a[], int ipiv[]) - { - int info; - SGETRF(&m,&m,a,&m,ipiv,&info); - for(int i = 0; i < m; ++i ){ - ipiv[i] -= 1; - } - } - - DLLEXPORT void d_lu_factor(int m, double a[], int ipiv[]) - { - int info; - DGETRF(&m,&m,a,&m,ipiv,&info); - for(int i = 0; i < m; ++i ){ - ipiv[i] -= 1; - } - } - - DLLEXPORT void c_lu_factor(int m, MKL_Complex8 a[], int ipiv[]) - { - int info; - CGETRF(&m,&m,a,&m,ipiv,&info); - for(int i = 0; i < m; ++i ){ - ipiv[i] -= 1; - } - } - - DLLEXPORT void z_lu_factor(int m, MKL_Complex16 a[], int ipiv[]) - { - int info; - ZGETRF(&m,&m,a,&m,ipiv,&info); - for(int i = 0; i < m; ++i ){ - ipiv[i] -= 1; - } - } - - DLLEXPORT void s_lu_inverse(int n, float a[], float work[], int lwork) - { - int* ipiv = new int[n]; - int info; - SGETRF(&n,&n,a,&n,ipiv,&info); - SGETRI(&n,a,&n,ipiv,work,&lwork,&info); + DLLEXPORT int s_cholesky_factor(int n, float a[]){ + char uplo = 'L'; + int info = 0; + SPOTRF(&uplo, &n, a, &n, &info); + for (int i = 0; i < n; ++i) + { + int index = i * n; + for (int j = 0; j < n && i > j; ++j) + { + a[index + j] = 0; + } + } + return info; + } + + DLLEXPORT int d_cholesky_factor(int n, double* a){ + char uplo = 'L'; + int info = 0; + DPOTRF(&uplo, &n, a, &n, &info); + for (int i = 0; i < n; ++i) + { + int index = i * n; + for (int j = 0; j < n && i > j; ++j) + { + a[index + j] = 0; + } + } + return info; + } + + DLLEXPORT int c_cholesky_factor(int n, Complex8 a[]){ + char uplo = 'L'; + int info = 0; + Complex8 zero; + zero.real = 0.0; + zero.real = 0.0; + CPOTRF(&uplo, &n, a, &n, &info); + for (int i = 0; i < n; ++i) + { + int index = i * n; + for (int j = 0; j < n && i > j; ++j) + { + a[index + j] = zero; + } + } + return info; + } + + DLLEXPORT int z_cholesky_factor(int n, Complex16 a[]){ + char uplo = 'L'; + int info = 0; + Complex16 zero; + zero.real = 0.0; + zero.real = 0.0; + ZPOTRF(&uplo, &n, a, &n, &info); + for (int i = 0; i < n; ++i) + { + int index = i * n; + for (int j = 0; j < n && i > j; ++j) + { + a[index + j] = zero; + } + } + return info; + } + + DLLEXPORT float s_matrix_norm(char norm, int m, int n, float a[], float work[]) + { + return SLANGE(&norm, &m, &n, a, &m, work); + } + + DLLEXPORT double d_matrix_norm(char norm, int m, int n, double a[], double work[]) + { + return DLANGE(&norm, &m, &n, a, &m, work); + } + + DLLEXPORT float c_matrix_norm(char norm, int m, int n, MKL_Complex8 a[], float work[]) + { + return CLANGE(&norm, &m, &n, a, &m, work); + } + + DLLEXPORT double z_matrix_norm(char norm, int m, int n, MKL_Complex16 a[], double work[]) + { + return ZLANGE(&norm, &m, &n, a, &m, work); + } + + DLLEXPORT void s_lu_factor(int m, float a[], int ipiv[]) + { + int info; + SGETRF(&m,&m,a,&m,ipiv,&info); + for(int i = 0; i < m; ++i ){ + ipiv[i] -= 1; + } + } + + DLLEXPORT void d_lu_factor(int m, double a[], int ipiv[]) + { + int info; + DGETRF(&m,&m,a,&m,ipiv,&info); + for(int i = 0; i < m; ++i ){ + ipiv[i] -= 1; + } + } + + DLLEXPORT void c_lu_factor(int m, MKL_Complex8 a[], int ipiv[]) + { + int info; + CGETRF(&m,&m,a,&m,ipiv,&info); + for(int i = 0; i < m; ++i ){ + ipiv[i] -= 1; + } + } + + DLLEXPORT void z_lu_factor(int m, MKL_Complex16 a[], int ipiv[]) + { + int info; + ZGETRF(&m,&m,a,&m,ipiv,&info); + for(int i = 0; i < m; ++i ){ + ipiv[i] -= 1; + } + } + + DLLEXPORT void s_lu_inverse(int n, float a[], float work[], int lwork) + { + int* ipiv = new int[n]; + int info; + SGETRF(&n,&n,a,&n,ipiv,&info); + SGETRI(&n,a,&n,ipiv,work,&lwork,&info); delete[] ipiv; - } - - DLLEXPORT void d_lu_inverse(int n, double a[], double work[], int lwork) - { - int* ipiv = new int[n]; - int info; - DGETRF(&n,&n,a,&n,ipiv,&info); - DGETRI(&n,a,&n,ipiv,work,&lwork,&info); + } + + DLLEXPORT void d_lu_inverse(int n, double a[], double work[], int lwork) + { + int* ipiv = new int[n]; + int info; + DGETRF(&n,&n,a,&n,ipiv,&info); + DGETRI(&n,a,&n,ipiv,work,&lwork,&info); delete[] ipiv; - } - - DLLEXPORT void c_lu_inverse(int n, MKL_Complex8 a[], MKL_Complex8 work[], int lwork) - { - int* ipiv = new int[n]; - int info; - CGETRF(&n,&n,a,&n,ipiv,&info); - CGETRI(&n,a,&n,ipiv,work,&lwork,&info); + } + + DLLEXPORT void c_lu_inverse(int n, MKL_Complex8 a[], MKL_Complex8 work[], int lwork) + { + int* ipiv = new int[n]; + int info; + CGETRF(&n,&n,a,&n,ipiv,&info); + CGETRI(&n,a,&n,ipiv,work,&lwork,&info); delete[] ipiv; - } - - DLLEXPORT void z_lu_inverse(int n, MKL_Complex16 a[], MKL_Complex16 work[], int lwork) - { - int* ipiv = new int[n]; - int info; - ZGETRF(&n,&n,a,&n,ipiv,&info); - ZGETRI(&n,a,&n,ipiv,work,&lwork,&info); + } + + DLLEXPORT void z_lu_inverse(int n, MKL_Complex16 a[], MKL_Complex16 work[], int lwork) + { + int* ipiv = new int[n]; + int info; + ZGETRF(&n,&n,a,&n,ipiv,&info); + ZGETRI(&n,a,&n,ipiv,work,&lwork,&info); delete[] ipiv; - } + } DLLEXPORT void s_lu_inverse_factored(int n, float a[], int ipiv[], float work[], int lwork) - { - int i; - for(i = 0; i < n; ++i ){ - ipiv[i] += 1; - } - int info; - SGETRI(&n,a,&n,ipiv,work,&lwork,&info); - - for(i = 0; i < n; ++i ){ - ipiv[i] -= 1; - } - } - - DLLEXPORT void d_lu_inverse_factored(int n, double a[], int ipiv[], double work[], int lwork) - { - int i; - for(i = 0; i < n; ++i ){ - ipiv[i] += 1; - } - - int info; - DGETRI(&n,a,&n,ipiv,work,&lwork,&info); - - for(i = 0; i < n; ++i ){ - ipiv[i] -= 1; - } - } - - DLLEXPORT void c_lu_inverse_factored(int n, MKL_Complex8 a[], int ipiv[], MKL_Complex8 work[], int lwork) - { - int i; - for(i = 0; i < n; ++i ){ - ipiv[i] += 1; - } - - int info; - CGETRI(&n,a,&n,ipiv,work,&lwork,&info); - - for(i = 0; i < n; ++i ){ - ipiv[i] -= 1; - } - } - - DLLEXPORT void z_lu_inverse_factored(int n, MKL_Complex16 a[], int ipiv[], MKL_Complex16 work[], int lwork) - { - int i; - for(i = 0; i < n; ++i ){ - ipiv[i] += 1; - } - - int info; - ZGETRI(&n,a,&n,ipiv,work,&lwork,&info); - - for(i = 0; i < n; ++i ){ - ipiv[i] -= 1; - } - } - - DLLEXPORT void s_lu_solve_factored(int n, int nrhs, float a[], int ipiv[], float b[]) - { - int info; - int i; - for(i = 0; i < n; ++i ){ - ipiv[i] += 1; - } - - char trans = 'N'; - SGETRS(&trans, &n, &nrhs, a, &n, ipiv, b, &n, &info); - for(i = 0; i < n; ++i ){ - ipiv[i] -= 1; - } - } - - DLLEXPORT void d_lu_solve_factored(int n, int nrhs, double a[], int ipiv[], double b[]) - { - int info; - int i; - for(i = 0; i < n; ++i ){ - ipiv[i] += 1; - } - - char trans = 'N'; - DGETRS(&trans, &n, &nrhs, a, &n, ipiv, b, &n, &info); - for(i = 0; i < n; ++i ){ - ipiv[i] -= 1; - } - } - - DLLEXPORT void c_lu_solve_factored(int n, int nrhs, MKL_Complex8 a[], int ipiv[], MKL_Complex8 b[]) - { - int info; - int i; - for(i = 0; i < n; ++i ){ - ipiv[i] += 1; - } - - char trans = 'N'; - CGETRS(&trans, &n, &nrhs, a, &n, ipiv, b, &n, &info); - for(i = 0; i < n; ++i ){ - ipiv[i] -= 1; - } - } - - DLLEXPORT void z_lu_solve_factored(int n, int nrhs, MKL_Complex16 a[], int ipiv[], MKL_Complex16 b[]) - { - int info; - int i; - for(i = 0; i < n; ++i ){ - ipiv[i] += 1; - } - - char trans = 'N'; - ZGETRS(&trans, &n, &nrhs, a, &n, ipiv, b, &n, &info); - for(i = 0; i < n; ++i ){ - ipiv[i] -= 1; - } - } + { + int i; + for(i = 0; i < n; ++i ){ + ipiv[i] += 1; + } + int info; + SGETRI(&n,a,&n,ipiv,work,&lwork,&info); + + for(i = 0; i < n; ++i ){ + ipiv[i] -= 1; + } + } + + DLLEXPORT void d_lu_inverse_factored(int n, double a[], int ipiv[], double work[], int lwork) + { + int i; + for(i = 0; i < n; ++i ){ + ipiv[i] += 1; + } + + int info; + DGETRI(&n,a,&n,ipiv,work,&lwork,&info); + + for(i = 0; i < n; ++i ){ + ipiv[i] -= 1; + } + } + + DLLEXPORT void c_lu_inverse_factored(int n, MKL_Complex8 a[], int ipiv[], MKL_Complex8 work[], int lwork) + { + int i; + for(i = 0; i < n; ++i ){ + ipiv[i] += 1; + } + + int info; + CGETRI(&n,a,&n,ipiv,work,&lwork,&info); + + for(i = 0; i < n; ++i ){ + ipiv[i] -= 1; + } + } + + DLLEXPORT void z_lu_inverse_factored(int n, MKL_Complex16 a[], int ipiv[], MKL_Complex16 work[], int lwork) + { + int i; + for(i = 0; i < n; ++i ){ + ipiv[i] += 1; + } + + int info; + ZGETRI(&n,a,&n,ipiv,work,&lwork,&info); + + for(i = 0; i < n; ++i ){ + ipiv[i] -= 1; + } + } + + DLLEXPORT void s_lu_solve_factored(int n, int nrhs, float a[], int ipiv[], float b[]) + { + int info; + int i; + for(i = 0; i < n; ++i ){ + ipiv[i] += 1; + } + + char trans ='N'; + SGETRS(&trans, &n, &nrhs, a, &n, ipiv, b, &n, &info); + for(i = 0; i < n; ++i ){ + ipiv[i] -= 1; + } + } + + DLLEXPORT void d_lu_solve_factored(int n, int nrhs, double a[], int ipiv[], double b[]) + { + int info; + int i; + for(i = 0; i < n; ++i ){ + ipiv[i] += 1; + } + + char trans ='N'; + DGETRS(&trans, &n, &nrhs, a, &n, ipiv, b, &n, &info); + for(i = 0; i < n; ++i ){ + ipiv[i] -= 1; + } + } + + DLLEXPORT void c_lu_solve_factored(int n, int nrhs, MKL_Complex8 a[], int ipiv[], MKL_Complex8 b[]) + { + int info; + int i; + for(i = 0; i < n; ++i ){ + ipiv[i] += 1; + } + + char trans ='N'; + CGETRS(&trans, &n, &nrhs, a, &n, ipiv, b, &n, &info); + for(i = 0; i < n; ++i ){ + ipiv[i] -= 1; + } + } + + DLLEXPORT void z_lu_solve_factored(int n, int nrhs, MKL_Complex16 a[], int ipiv[], MKL_Complex16 b[]) + { + int info; + int i; + for(i = 0; i < n; ++i ){ + ipiv[i] += 1; + } + + char trans ='N'; + ZGETRS(&trans, &n, &nrhs, a, &n, ipiv, b, &n, &info); + for(i = 0; i < n; ++i ){ + ipiv[i] -= 1; + } + } DLLEXPORT void s_lu_solve(int n, int nrhs, float a[], float b[]) - { - int* ipiv = new int[n]; - int info; - SGETRF(&n,&n,a,&n,ipiv,&info); - - char trans = 'N'; - SGETRS(&trans, &n, &nrhs, a, &n, ipiv, b, &n, &info); + { + float* clone = new float[n*n]; + std::memcpy(clone, a, n*n*sizeof(float)); + + int* ipiv = new int[n]; + int info; + SGETRF(&n, &n, clone, &n, ipiv, &info); + char trans ='N'; + SGETRS(&trans, &n, &nrhs, clone, &n, ipiv, b, &n, &info); delete[] ipiv; - } - - DLLEXPORT void d_lu_solve(int n, int nrhs, double a[], double b[]) - { - int* ipiv = new int[n]; - int info; - DGETRF(&n,&n,a,&n,ipiv,&info); - - char trans = 'N'; - DGETRS(&trans, &n, &nrhs, a, &n, ipiv, b, &n, &info); + delete[] clone; + } + + DLLEXPORT void d_lu_solve(int n, int nrhs, double a[], double b[]) + { + double* clone = new double[n*n]; + std::memcpy(clone, a, n*n*sizeof(double)); + + int* ipiv = new int[n]; + int info; + DGETRF(&n, &n, clone, &n, ipiv, &info); + char trans ='N'; + DGETRS(&trans, &n, &nrhs, clone, &n, ipiv, b, &n, &info); delete[] ipiv; - } - - DLLEXPORT void c_lu_solve(int n, int nrhs, MKL_Complex8 a[], MKL_Complex8 b[]) - { - int* ipiv = new int[n]; - int info; - CGETRF(&n,&n,a,&n,ipiv,&info); - - char trans = 'N'; - CGETRS(&trans, &n, &nrhs, a, &n, ipiv, b, &n, &info); + delete[] clone; + } + + DLLEXPORT void c_lu_solve(int n, int nrhs, MKL_Complex8 a[], MKL_Complex8 b[]) + { + MKL_Complex8* clone = new MKL_Complex8[n*n]; + std::memcpy(clone, a, n*n*sizeof(MKL_Complex8)); + + int* ipiv = new int[n]; + int info; + CGETRF(&n, &n, clone, &n, ipiv, &info); + char trans ='N'; + CGETRS(&trans, &n, &nrhs, clone, &n, ipiv, b, &n, &info); delete[] ipiv; - } - - DLLEXPORT void z_lu_solve(int n, int nrhs, MKL_Complex16 a[], MKL_Complex16 b[]) - { - int* ipiv = new int[n]; - int info; - ZGETRF(&n,&n,a,&n,ipiv,&info); - - char trans = 'N'; - ZGETRS(&trans, &n, &nrhs, a, &n, ipiv, b, &n, &info); + delete[] clone; + } + + DLLEXPORT void z_lu_solve(int n, int nrhs, MKL_Complex16 a[], MKL_Complex16 b[]) + { + MKL_Complex16* clone = new MKL_Complex16[n*n]; + std::memcpy(clone, a, n*n*sizeof(MKL_Complex16)); + + int* ipiv = new int[n]; + int info; + ZGETRF(&n, &n, clone, &n, ipiv, &info); + char trans ='N'; + ZGETRS(&trans, &n, &nrhs, clone, &n, ipiv, b, &n, &info); delete[] ipiv; - } + delete[] clone; + } - DLLEXPORT void s_cholesky_solve(int n, int nrhs, float a[], float b[]) - { - char uplo = 'L'; - int info = 0; - SPOTRS(&uplo, &n, &nrhs, a, &n, b, &n, &info); - } - - DLLEXPORT void d_cholesky_solve(int n, int nrhs, double a[], double b[]) - { - char uplo = 'L'; - int info = 0; - DPOTRS(&uplo, &n, &nrhs, a, &n, b, &n, &info); - } - - DLLEXPORT void c_cholesky_solve(int n, int nrhs, MKL_Complex8 a[], MKL_Complex8 b[]) - { - char uplo = 'L'; - int info = 0; - CPOTRS(&uplo, &n, &nrhs, a, &n, b, &n, &info); - } - - DLLEXPORT void z_cholesky_solve(int n, int nrhs, MKL_Complex16 a[], MKL_Complex16 b[]) - { - char uplo = 'L'; - int info = 0; - ZPOTRS(&uplo, &n, &nrhs, a, &n, b, &n, &info); - } - - DLLEXPORT void s_qr_factor(int m, int n, float r[], float tau[], float q[], float work[], int len) - { - int info = 0; - SGEQRF(&m, &n, r, &m, tau, work, &len, &info); - - for (int i = 0; i < m; ++i) - { - for (int j = 0; j < m && j < n; ++j) - { - if (i > j) - { - q[j * m + i] = r[j * m + i]; - } - } - } - - //compute the q elements explicitly - if (m <= n) - { - SORGQR(&m, &m, &m, q, &m, tau, work, &len, &info); - } - else - { - SORGQR(&m, &n, &n, q, &m, tau, work, &len, &info); - } - } - - DLLEXPORT void d_qr_factor(int m, int n, double r[], double tau[], double q[], double work[], int len) - { - int info = 0; - DGEQRF(&m, &n, r, &m, tau, work, &len, &info); - - for (int i = 0; i < m; ++i) - { - for (int j = 0; j < m && j < n; ++j) - { - if (i > j) - { - q[j * m + i] = r[j * m + i]; - } - } - } - - //compute the q elements explicitly - if (m <= n) - { - DORGQR(&m, &m, &m, q, &m, tau, work, &len, &info); - } - else - { - DORGQR(&m, &n, &n, q, &m, tau, work, &len, &info); - } - } - - DLLEXPORT void c_qr_factor(int m, int n, MKL_Complex8 r[], MKL_Complex8 tau[], MKL_Complex8 q[], MKL_Complex8 work[], int len) - { - int info = 0; - CGEQRF(&m, &n, r, &m, tau, work, &len, &info); - - for (int i = 0; i < m; ++i) - { - for (int j = 0; j < m && j < n; ++j) - { - if (i > j) - { - q[j * m + i] = r[j * m + i]; - } - } - } - - //compute the q elements explicitly - if (m <= n) - { - CUNGQR(&m, &m, &m, q, &m, tau, work, &len, &info); - } - else - { - CUNGQR(&m, &n, &n, q, &m, tau, work, &len, &info); - } - } - - DLLEXPORT void z_qr_factor(int m, int n, MKL_Complex16 r[], MKL_Complex16 tau[], MKL_Complex16 q[], MKL_Complex16 work[], int len) - { - int info = 0; - ZGEQRF(&m, &n, r, &m, tau, work, &len, &info); - - for (int i = 0; i < m; ++i) - { - for (int j = 0; j < m && j < n; ++j) - { - if (i > j) - { - q[j * m + i] = r[j * m + i]; - } - } - } - - //compute the q elements explicitly - if (m <= n) - { - ZUNGQR(&m, &m, &m, q, &m, tau, work, &len, &info); - } - else - { - ZUNGQR(&m, &n, &n, q, &m, tau, work, &len, &info); - } - } - - DLLEXPORT void s_qr_solve(int m, int n, int bn, float r[], float b[], float tau[], float x[], float work[], int len) - { - char side ='L'; - char tran = 'T'; - int info = 0; - SORMQR(&side, &tran, &m, &bn, &n, r, &m, tau, b, &m, work, &len, &info); - cblas_strsm(CblasColMajor,CblasLeft,CblasUpper,CblasNoTrans,CblasNonUnit, n, bn, 1.0, r, m, b, m); - for (int i = 0; i < n; ++i) - { - for (int j = 0; j < bn; ++j) - { - x[j * n + i] = b[j * m + i]; - } - } - } - - DLLEXPORT void d_qr_solve(int m, int n, int bn, double r[], double b[], double tau[], double x[], double work[], int len) - { - char side ='L'; - char tran = 'T'; - int info = 0; - DORMQR(&side, &tran, &m, &bn, &n, r, &m, tau, b, &m, work, &len, &info); - cblas_dtrsm(CblasColMajor,CblasLeft,CblasUpper,CblasNoTrans,CblasNonUnit, n, bn, 1.0, r, m, b, m); - for (int i = 0; i < n; ++i) - { - for (int j = 0; j < bn; ++j) - { - x[j * n + i] = b[j * m + i]; - } - } - } - - DLLEXPORT void c_qr_solve(int m, int n, int bn, MKL_Complex8 r[], MKL_Complex8 b[], MKL_Complex8 tau[], MKL_Complex8 x[], MKL_Complex8 work[], int len) - { - char side ='L'; - char tran = 'T'; - int info = 0; - CUNMQR(&side, &tran, &m, &bn, &n, r, &m, tau, b, &m, work, &len, &info); - MKL_Complex8 one; - one.real = 1.0; - cblas_ctrsm(CblasColMajor,CblasLeft,CblasUpper,CblasNoTrans,CblasNonUnit, n, bn, &one, r, m, b, m); - for (int i = 0; i < n; ++i) - { - for (int j = 0; j < bn; ++j) - { - x[j * n + i] = b[j * m + i]; - } - } - } - - DLLEXPORT void z_qr_solve(int m, int n, int bn, MKL_Complex16 r[], MKL_Complex16 b[], MKL_Complex16 tau[], MKL_Complex16 x[], MKL_Complex16 work[], int len) - { - char side ='L'; - char tran = 'T'; - int info = 0; - ZUNMQR(&side, &tran, &m, &bn, &n, r, &m, tau, b, &m, work, &len, &info); - MKL_Complex16 one; - one.real = 1.0; - cblas_ztrsm(CblasColMajor,CblasLeft,CblasUpper,CblasNoTrans,CblasNonUnit, n, bn, &one, r, m, b, m); - for (int i = 0; i < n; ++i) - { - for (int j = 0; j < bn; ++j) - { - x[j * n + i] = b[j * m + i]; - } - } - } - - DLLEXPORT void s_svd_factor(bool compute_vectors, int m, int n, float a[], float s[], float u[], float v[], float work[], int len) - { - int info = 0; - char job = compute_vectors ? 'A' : 'N'; - SGESVD(&job, &job, &m, &n, a, &m, s, u, &m, v, &n, work, &len, &info); - } - - DLLEXPORT void d_svd_factor(bool compute_vectors, int m, int n, double a[], double s[], double u[], double v[], double work[], int len) - { - int info = 0; - char job = compute_vectors ? 'A' : 'N'; - DGESVD(&job, &job, &m, &n, a, &m, s, u, &m, v, &n, work, &len, &info); - } - - DLLEXPORT void c_svd_factor(bool compute_vectors, int m, int n, MKL_Complex8 a[], float s[], MKL_Complex8 u[], MKL_Complex8 v[], MKL_Complex8 work[], int len, float rwork[]) - { - int info = 0; - char job = compute_vectors ? 'A' : 'N'; - CGESVD(&job, &job, &m, &n, a, &m, s, u, &m, v, &n, work, &len, rwork, &info); - } - - DLLEXPORT void z_svd_factor(bool compute_vectors, int m, int n, MKL_Complex16 a[], double s[], MKL_Complex16 u[], MKL_Complex16 v[], MKL_Complex16 work[], int len, double rwork[]) - { - int info = 0; - char job = compute_vectors ? 'A' : 'N'; - ZGESVD(&job, &job, &m, &n, a, &m, s, u, &m, v, &n, work, &len, rwork, &info); - } + { + char uplo = 'L'; + int info = 0; + SPOTRS(&uplo, &n, &nrhs, a, &n, b, &n, &info); + } + + DLLEXPORT void d_cholesky_solve(int n, int nrhs, double a[], double b[]) + { + char uplo = 'L'; + int info = 0; + DPOTRS(&uplo, &n, &nrhs, a, &n, b, &n, &info); + } + + DLLEXPORT void c_cholesky_solve(int n, int nrhs, MKL_Complex8 a[], MKL_Complex8 b[]) + { + char uplo = 'L'; + int info = 0; + CPOTRS(&uplo, &n, &nrhs, a, &n, b, &n, &info); + } + + DLLEXPORT void z_cholesky_solve(int n, int nrhs, MKL_Complex16 a[], MKL_Complex16 b[]) + { + char uplo = 'L'; + int info = 0; + ZPOTRS(&uplo, &n, &nrhs, a, &n, b, &n, &info); + } + + DLLEXPORT void s_qr_factor(int m, int n, float r[], float tau[], float q[], float work[], int len) + { + int info = 0; + SGEQRF(&m, &n, r, &m, tau, work, &len, &info); + + for (int i = 0; i < m; ++i) + { + for (int j = 0; j < m && j < n; ++j) + { + if (i > j) + { + q[j * m + i] = r[j * m + i]; + } + } + } + + //compute the q elements explicitly + if (m <= n) + { + SORGQR(&m, &m, &m, q, &m, tau, work, &len, &info); + } + else + { + SORGQR(&m, &n, &n, q, &m, tau, work, &len, &info); + } + } + + DLLEXPORT void d_qr_factor(int m, int n, double r[], double tau[], double q[], double work[], int len) + { + int info = 0; + DGEQRF(&m, &n, r, &m, tau, work, &len, &info); + + for (int i = 0; i < m; ++i) + { + for (int j = 0; j < m && j < n; ++j) + { + if (i > j) + { + q[j * m + i] = r[j * m + i]; + } + } + } + + //compute the q elements explicitly + if (m <= n) + { + DORGQR(&m, &m, &m, q, &m, tau, work, &len, &info); + } + else + { + DORGQR(&m, &n, &n, q, &m, tau, work, &len, &info); + } + } + + DLLEXPORT void c_qr_factor(int m, int n, MKL_Complex8 r[], MKL_Complex8 tau[], MKL_Complex8 q[], MKL_Complex8 work[], int len) + { + int info = 0; + CGEQRF(&m, &n, r, &m, tau, work, &len, &info); + + for (int i = 0; i < m; ++i) + { + for (int j = 0; j < m && j < n; ++j) + { + if (i > j) + { + q[j * m + i] = r[j * m + i]; + } + } + } + + //compute the q elements explicitly + if (m <= n) + { + CUNGQR(&m, &m, &m, q, &m, tau, work, &len, &info); + } + else + { + CUNGQR(&m, &n, &n, q, &m, tau, work, &len, &info); + } + } + + DLLEXPORT void z_qr_factor(int m, int n, MKL_Complex16 r[], MKL_Complex16 tau[], MKL_Complex16 q[], MKL_Complex16 work[], int len) + { + int info = 0; + ZGEQRF(&m, &n, r, &m, tau, work, &len, &info); + + for (int i = 0; i < m; ++i) + { + for (int j = 0; j < m && j < n; ++j) + { + if (i > j) + { + q[j * m + i] = r[j * m + i]; + } + } + } + + //compute the q elements explicitly + if (m <= n) + { + ZUNGQR(&m, &m, &m, q, &m, tau, work, &len, &info); + } + else + { + ZUNGQR(&m, &n, &n, q, &m, tau, work, &len, &info); + } + } + + DLLEXPORT void s_qr_solve(int m, int n, int bn, float r[], float b[], float tau[], float x[], float work[], int len) + { + char side ='L'; + char tran = 'T'; + int info = 0; + SORMQR(&side, &tran, &m, &bn, &n, r, &m, tau, b, &m, work, &len, &info); + cblas_strsm(CblasColMajor,CblasLeft,CblasUpper,CblasNoTrans,CblasNonUnit, n, bn, 1.0, r, m, b, m); + for (int i = 0; i < n; ++i) + { + for (int j = 0; j < bn; ++j) + { + x[j * n + i] = b[j * m + i]; + } + } + } + + DLLEXPORT void d_qr_solve(int m, int n, int bn, double r[], double b[], double tau[], double x[], double work[], int len) + { + char side ='L'; + char tran = 'T'; + int info = 0; + DORMQR(&side, &tran, &m, &bn, &n, r, &m, tau, b, &m, work, &len, &info); + cblas_dtrsm(CblasColMajor,CblasLeft,CblasUpper,CblasNoTrans,CblasNonUnit, n, bn, 1.0, r, m, b, m); + for (int i = 0; i < n; ++i) + { + for (int j = 0; j < bn; ++j) + { + x[j * n + i] = b[j * m + i]; + } + } + } + + DLLEXPORT void c_qr_solve(int m, int n, int bn, MKL_Complex8 r[], MKL_Complex8 b[], MKL_Complex8 tau[], MKL_Complex8 x[], MKL_Complex8 work[], int len) + { + char side ='L'; + char tran = 'T'; + int info = 0; + CUNMQR(&side, &tran, &m, &bn, &n, r, &m, tau, b, &m, work, &len, &info); + MKL_Complex8 one; + one.real = 1.0; + cblas_ctrsm(CblasColMajor,CblasLeft,CblasUpper,CblasNoTrans,CblasNonUnit, n, bn, &one, r, m, b, m); + for (int i = 0; i < n; ++i) + { + for (int j = 0; j < bn; ++j) + { + x[j * n + i] = b[j * m + i]; + } + } + } + + DLLEXPORT void z_qr_solve(int m, int n, int bn, MKL_Complex16 r[], MKL_Complex16 b[], MKL_Complex16 tau[], MKL_Complex16 x[], MKL_Complex16 work[], int len) + { + char side ='L'; + char tran = 'T'; + int info = 0; + ZUNMQR(&side, &tran, &m, &bn, &n, r, &m, tau, b, &m, work, &len, &info); + MKL_Complex16 one; + one.real = 1.0; + cblas_ztrsm(CblasColMajor,CblasLeft,CblasUpper,CblasNoTrans,CblasNonUnit, n, bn, &one, r, m, b, m); + for (int i = 0; i < n; ++i) + { + for (int j = 0; j < bn; ++j) + { + x[j * n + i] = b[j * m + i]; + } + } + } + + DLLEXPORT void s_svd_factor(bool compute_vectors, int m, int n, float a[], float s[], float u[], float v[], float work[], int len) + { + int info = 0; + char job = compute_vectors ? 'A' : 'N'; + SGESVD(&job, &job, &m, &n, a, &m, s, u, &m, v, &n, work, &len, &info); + } + + DLLEXPORT void d_svd_factor(bool compute_vectors, int m, int n, double a[], double s[], double u[], double v[], double work[], int len) + { + int info = 0; + char job = compute_vectors ? 'A' : 'N'; + DGESVD(&job, &job, &m, &n, a, &m, s, u, &m, v, &n, work, &len, &info); + } + + DLLEXPORT void c_svd_factor(bool compute_vectors, int m, int n, MKL_Complex8 a[], float s[], MKL_Complex8 u[], MKL_Complex8 v[], MKL_Complex8 work[], int len, float rwork[]) + { + int info = 0; + char job = compute_vectors ? 'A' : 'N'; + CGESVD(&job, &job, &m, &n, a, &m, s, u, &m, v, &n, work, &len, rwork, &info); + } + + DLLEXPORT void z_svd_factor(bool compute_vectors, int m, int n, MKL_Complex16 a[], double s[], MKL_Complex16 u[], MKL_Complex16 v[], MKL_Complex16 work[], int len, double rwork[]) + { + int info = 0; + char job = compute_vectors ? 'A' : 'N'; + ZGESVD(&job, &job, &m, &n, a, &m, s, u, &m, v, &n, work, &len, rwork, &info); + } } \ No newline at end of file diff --git a/src/Numerics/Algorithms/LinearAlgebra/ILinearAlgebraProviderOfT.cs b/src/Numerics/Algorithms/LinearAlgebra/ILinearAlgebraProviderOfT.cs index 43de7f26..08bea682 100644 --- a/src/Numerics/Algorithms/LinearAlgebra/ILinearAlgebraProviderOfT.cs +++ b/src/Numerics/Algorithms/LinearAlgebra/ILinearAlgebraProviderOfT.cs @@ -3,7 +3,7 @@ // http://numerics.mathdotnet.com // http://github.com/mathnet/mathnet-numerics // http://mathnetnumerics.codeplex.com -// Copyright (c) 2009-2010 Math.NET +// Copyright (c) 2009-2011 Math.NET // Permission is hereby granted, free of charge, to any person // obtaining a copy of this software and associated documentation // files (the "Software"), to deal in the Software without @@ -283,29 +283,6 @@ namespace MathNet.Numerics.Algorithms.LinearAlgebra /// This is equivalent to the GETRS LAPACK routine. void LUSolveFactored(int columnsOfB, T[] a, int order, int[] ipiv, T[] b); - /// - /// Solves A*X=B for X using LU factorization. - /// - /// How to transpose the matrix. - /// The number of columns of B. - /// The square matrix A. - /// The order of the square matrix . - /// The B matrix. - /// This is equivalent to the GETRF and GETRS LAPACK routines. - void LUSolve(Transpose transposeA, int columnsOfB, T[] a, int order, T[] b); - - /// - /// Solves A*X=B for X using a previously factored A matrix. - /// - /// How to transpose the matrix. - /// The number of columns of B. - /// The factored A matrix. - /// The order of the square matrix . - /// The pivot indices of . - /// The B matrix. - /// This is equivalent to the GETRS LAPACK routine. - void LUSolveFactored(Transpose transposeA, int columnsOfB, T[] a, int order, int[] ipiv, T[] b); - /// /// Computes the Cholesky factorization of A. /// diff --git a/src/Numerics/Algorithms/LinearAlgebra/ManagedLinearAlgebraProvider.Complex.cs b/src/Numerics/Algorithms/LinearAlgebra/ManagedLinearAlgebraProvider.Complex.cs index a8246018..43027bba 100644 --- a/src/Numerics/Algorithms/LinearAlgebra/ManagedLinearAlgebraProvider.Complex.cs +++ b/src/Numerics/Algorithms/LinearAlgebra/ManagedLinearAlgebraProvider.Complex.cs @@ -1026,7 +1026,7 @@ namespace MathNet.Numerics.Algorithms.LinearAlgebra /// The number of columns of B. /// The square matrix A. /// The order of the square matrix . - /// The B matrix. + /// On input the B matrix; on output the X matrix. /// This is equivalent to the GETRF and GETRS LAPACK routines. public virtual void LUSolve(int columnsOfB, Complex[] a, int order, Complex[] b) { @@ -1051,8 +1051,10 @@ namespace MathNet.Numerics.Algorithms.LinearAlgebra } var ipiv = new int[order]; - LUFactor(a, order, ipiv); - LUSolveFactored(columnsOfB, a, order, ipiv, b); + var clone = new Complex[a.Length]; + Array.Copy(a, 0, clone, 0, a.Length); + LUFactor(clone, order, ipiv); + LUSolveFactored(columnsOfB, clone, order, ipiv, b); } /// @@ -1062,7 +1064,7 @@ namespace MathNet.Numerics.Algorithms.LinearAlgebra /// The factored A matrix. /// The order of the square matrix . /// The pivot indices of . - /// The B matrix. + /// On input the B matrix; on output the X matrix. /// This is equivalent to the GETRS LAPACK routine. public virtual void LUSolveFactored(int columnsOfB, Complex[] a, int order, int[] ipiv, Complex[] b) { @@ -1151,116 +1153,6 @@ namespace MathNet.Numerics.Algorithms.LinearAlgebra } } - /// - /// Solves A*X=B for X using LU factorization. - /// - /// How to transpose the matrix. - /// The number of columns of B. - /// The square matrix A. - /// The order of the square matrix . - /// The B matrix. - /// This is equivalent to the GETRF and GETRS LAPACK routines. - public virtual void LUSolve(Transpose transposeA, int columnsOfB, Complex[] a, int order, Complex[] b) - { - if (a == null) - { - throw new ArgumentNullException("a"); - } - - if (b == null) - { - throw new ArgumentNullException("b"); - } - - if (a.Length != order * order) - { - throw new ArgumentException(Resources.ArgumentArraysSameLength, "a"); - } - - if (b.Length != order * columnsOfB) - { - throw new ArgumentException(Resources.ArgumentArraysSameLength, "b"); - } - - var ipiv = new int[order]; - LUFactor(a, order, ipiv); - LUSolveFactored(transposeA, columnsOfB, a, order, ipiv, b); - } - - /// - /// Solves A*X=B for X using a previously factored A matrix. - /// - /// How to transpose the matrix. - /// The number of columns of B. - /// The factored A matrix. - /// The order of the square matrix . - /// The pivot indices of . - /// The B matrix. - /// This is equivalent to the GETRS LAPACK routine. - public virtual void LUSolveFactored(Transpose transposeA, int columnsOfB, Complex[] a, int order, int[] ipiv, Complex[] b) - { - if (a == null) - { - throw new ArgumentNullException("a"); - } - - if (ipiv == null) - { - throw new ArgumentNullException("ipiv"); - } - - if (b == null) - { - throw new ArgumentNullException("b"); - } - - if (a.Length != order * order) - { - throw new ArgumentException(Resources.ArgumentArraysSameLength, "a"); - } - - if (ipiv.Length != order) - { - throw new ArgumentException(Resources.ArgumentArraysSameLength, "ipiv"); - } - - if (b.Length != order * columnsOfB) - { - throw new ArgumentException(Resources.ArgumentArraysSameLength, "b"); - } - - if (transposeA == Transpose.Transpose) - { - var aT = new Complex[a.Length]; - for (var i = 0; i < order; i++) - { - for (var j = 0; j < order; j++) - { - aT[(j * order) + i] = a[(i * order) + j]; - } - } - - LUSolveFactored(columnsOfB, aT, order, ipiv, b); - } - else if (transposeA == Transpose.ConjugateTranspose) - { - var acT = new Complex[a.Length]; - for (var i = 0; i < order; i++) - { - for (var j = 0; j < order; j++) - { - acT[(j * order) + i] = a[(i * order) + j].Conjugate(); - } - } - - LUSolveFactored(columnsOfB, acT, order, ipiv, b); - } - else - { - LUSolveFactored(columnsOfB, a, order, ipiv, b); - } - } - /// /// Computes the Cholesky factorization of A. /// diff --git a/src/Numerics/Algorithms/LinearAlgebra/ManagedLinearAlgebraProvider.Complex32.cs b/src/Numerics/Algorithms/LinearAlgebra/ManagedLinearAlgebraProvider.Complex32.cs index 40ba3cc6..e284ba47 100644 --- a/src/Numerics/Algorithms/LinearAlgebra/ManagedLinearAlgebraProvider.Complex32.cs +++ b/src/Numerics/Algorithms/LinearAlgebra/ManagedLinearAlgebraProvider.Complex32.cs @@ -1026,7 +1026,7 @@ namespace MathNet.Numerics.Algorithms.LinearAlgebra /// The number of columns of B. /// The square matrix A. /// The order of the square matrix . - /// The B matrix. + /// On input the B matrix; on output the X matrix. /// This is equivalent to the GETRF and GETRS LAPACK routines. public virtual void LUSolve(int columnsOfB, Complex32[] a, int order, Complex32[] b) { @@ -1051,8 +1051,10 @@ namespace MathNet.Numerics.Algorithms.LinearAlgebra } var ipiv = new int[order]; - LUFactor(a, order, ipiv); - LUSolveFactored(columnsOfB, a, order, ipiv, b); + var clone = new Complex32[a.Length]; + Array.Copy(a, 0, clone, 0, a.Length); + LUFactor(clone, order, ipiv); + LUSolveFactored(columnsOfB, clone, order, ipiv, b); } /// @@ -1062,7 +1064,7 @@ namespace MathNet.Numerics.Algorithms.LinearAlgebra /// The factored A matrix. /// The order of the square matrix . /// The pivot indices of . - /// The B matrix. + /// On input the B matrix; on output the X matrix. /// This is equivalent to the GETRS LAPACK routine. public virtual void LUSolveFactored(int columnsOfB, Complex32[] a, int order, int[] ipiv, Complex32[] b) { @@ -1151,116 +1153,6 @@ namespace MathNet.Numerics.Algorithms.LinearAlgebra } } - /// - /// Solves A*X=B for X using LU factorization. - /// - /// How to transpose the matrix. - /// The number of columns of B. - /// The square matrix A. - /// The order of the square matrix . - /// The B matrix. - /// This is equivalent to the GETRF and GETRS LAPACK routines. - public virtual void LUSolve(Transpose transposeA, int columnsOfB, Complex32[] a, int order, Complex32[] b) - { - if (a == null) - { - throw new ArgumentNullException("a"); - } - - if (b == null) - { - throw new ArgumentNullException("b"); - } - - if (a.Length != order * order) - { - throw new ArgumentException(Resources.ArgumentArraysSameLength, "a"); - } - - if (b.Length != order * columnsOfB) - { - throw new ArgumentException(Resources.ArgumentArraysSameLength, "b"); - } - - var ipiv = new int[order]; - LUFactor(a, order, ipiv); - LUSolveFactored(transposeA, columnsOfB, a, order, ipiv, b); - } - - /// - /// Solves A*X=B for X using a previously factored A matrix. - /// - /// How to transpose the matrix. - /// The number of columns of B. - /// The factored A matrix. - /// The order of the square matrix . - /// The pivot indices of . - /// The B matrix. - /// This is equivalent to the GETRS LAPACK routine. - public virtual void LUSolveFactored(Transpose transposeA, int columnsOfB, Complex32[] a, int order, int[] ipiv, Complex32[] b) - { - if (a == null) - { - throw new ArgumentNullException("a"); - } - - if (ipiv == null) - { - throw new ArgumentNullException("ipiv"); - } - - if (b == null) - { - throw new ArgumentNullException("b"); - } - - if (a.Length != order * order) - { - throw new ArgumentException(Resources.ArgumentArraysSameLength, "a"); - } - - if (ipiv.Length != order) - { - throw new ArgumentException(Resources.ArgumentArraysSameLength, "ipiv"); - } - - if (b.Length != order * columnsOfB) - { - throw new ArgumentException(Resources.ArgumentArraysSameLength, "b"); - } - - if (transposeA == Transpose.Transpose) - { - var aT = new Complex32[a.Length]; - for (var i = 0; i < order; i++) - { - for (var j = 0; j < order; j++) - { - aT[(j * order) + i] = a[(i * order) + j]; - } - } - - LUSolveFactored(columnsOfB, aT, order, ipiv, b); - } - else if (transposeA == Transpose.ConjugateTranspose) - { - var acT = new Complex32[a.Length]; - for (var i = 0; i < order; i++) - { - for (var j = 0; j < order; j++) - { - acT[(j * order) + i] = a[(i * order) + j].Conjugate(); - } - } - - LUSolveFactored(columnsOfB, acT, order, ipiv, b); - } - else - { - LUSolveFactored(columnsOfB, a, order, ipiv, b); - } - } - /// /// Computes the Cholesky factorization of A. /// diff --git a/src/Numerics/Algorithms/LinearAlgebra/ManagedLinearAlgebraProvider.Double.cs b/src/Numerics/Algorithms/LinearAlgebra/ManagedLinearAlgebraProvider.Double.cs index 4201a22d..ea906a31 100644 --- a/src/Numerics/Algorithms/LinearAlgebra/ManagedLinearAlgebraProvider.Double.cs +++ b/src/Numerics/Algorithms/LinearAlgebra/ManagedLinearAlgebraProvider.Double.cs @@ -1020,7 +1020,7 @@ namespace MathNet.Numerics.Algorithms.LinearAlgebra /// The number of columns of B. /// The square matrix A. /// The order of the square matrix . - /// The B matrix. + /// On input the B matrix; on output the X matrix. /// This is equivalent to the GETRF and GETRS LAPACK routines. public virtual void LUSolve(int columnsOfB, double[] a, int order, double[] b) { @@ -1045,8 +1045,10 @@ namespace MathNet.Numerics.Algorithms.LinearAlgebra } var ipiv = new int[order]; - LUFactor(a, order, ipiv); - LUSolveFactored(columnsOfB, a, order, ipiv, b); + var clone = new double[a.Length]; + Buffer.BlockCopy(a, 0, clone, 0, a.Length * Constants.SizeOfDouble); + LUFactor(clone, order, ipiv); + LUSolveFactored(columnsOfB, clone, order, ipiv, b); } /// @@ -1056,7 +1058,7 @@ namespace MathNet.Numerics.Algorithms.LinearAlgebra /// The factored A matrix. /// The order of the square matrix . /// The pivot indices of . - /// The B matrix. + /// On input the B matrix; on output the X matrix. /// This is equivalent to the GETRS LAPACK routine. public virtual void LUSolveFactored(int columnsOfB, double[] a, int order, int[] ipiv, double[] b) { @@ -1145,103 +1147,6 @@ namespace MathNet.Numerics.Algorithms.LinearAlgebra } } - /// - /// Solves A*X=B for X using LU factorization. - /// - /// How to transpose the matrix. - /// The number of columns of B. - /// The square matrix A. - /// The order of the square matrix . - /// The B matrix. - /// This is equivalent to the GETRF and GETRS LAPACK routines. - public virtual void LUSolve(Transpose transposeA, int columnsOfB, double[] a, int order, double[] b) - { - if (a == null) - { - throw new ArgumentNullException("a"); - } - - if (b == null) - { - throw new ArgumentNullException("b"); - } - - if (a.Length != order * order) - { - throw new ArgumentException(Resources.ArgumentArraysSameLength, "a"); - } - - if (b.Length != order * columnsOfB) - { - throw new ArgumentException(Resources.ArgumentArraysSameLength, "b"); - } - - var ipiv = new int[order]; - LUFactor(a, order, ipiv); - LUSolveFactored(transposeA, columnsOfB, a, order, ipiv, b); - } - - /// - /// Solves A*X=B for X using a previously factored A matrix. - /// - /// How to transpose the matrix. - /// The number of columns of B. - /// The factored A matrix. - /// The order of the square matrix . - /// The pivot indices of . - /// The B matrix. - /// This is equivalent to the GETRS LAPACK routine. - public virtual void LUSolveFactored(Transpose transposeA, int columnsOfB, double[] a, int order, int[] ipiv, double[] b) - { - if (a == null) - { - throw new ArgumentNullException("a"); - } - - if (ipiv == null) - { - throw new ArgumentNullException("ipiv"); - } - - if (b == null) - { - throw new ArgumentNullException("b"); - } - - if (a.Length != order * order) - { - throw new ArgumentException(Resources.ArgumentArraysSameLength, "a"); - } - - if (ipiv.Length != order) - { - throw new ArgumentException(Resources.ArgumentArraysSameLength, "ipiv"); - } - - if (b.Length != order * columnsOfB) - { - throw new ArgumentException(Resources.ArgumentArraysSameLength, "b"); - } - - if ((transposeA == Transpose.Transpose) || (transposeA == Transpose.ConjugateTranspose)) - { - var aT = new double[a.Length]; - for (var i = 0; i < order; i++) - { - for (var j = 0; j < order; j++) - { - aT[(j * order) + i] = a[(i * order) + j]; - } - } - - LUSolveFactored(columnsOfB, aT, order, ipiv, b); - } - else - { - LUSolveFactored(columnsOfB, a, order, ipiv, b); - } - } - /// /// Computes the Cholesky factorization of A. /// @@ -1829,7 +1734,7 @@ namespace MathNet.Numerics.Algorithms.LinearAlgebra var im = i * rowsR; var sum = 0.0; - for (var k = 0; k < rowsR; k++ ) + for (var k = 0; k < rowsR; k++) { sum += q[im + k] * column[k]; } @@ -2021,7 +1926,6 @@ namespace MathNet.Numerics.Algorithms.LinearAlgebra { // Compute the transformation for the l-th column and // place the l-th diagonal in vector s[l]. - var sum = 0.0; for (var i1 = l; i1 < rowsA; i1++) { diff --git a/src/Numerics/Algorithms/LinearAlgebra/ManagedLinearAlgebraProvider.Single.cs b/src/Numerics/Algorithms/LinearAlgebra/ManagedLinearAlgebraProvider.Single.cs index 9358cc19..8a125cc7 100644 --- a/src/Numerics/Algorithms/LinearAlgebra/ManagedLinearAlgebraProvider.Single.cs +++ b/src/Numerics/Algorithms/LinearAlgebra/ManagedLinearAlgebraProvider.Single.cs @@ -1021,7 +1021,7 @@ namespace MathNet.Numerics.Algorithms.LinearAlgebra /// The number of columns of B. /// The square matrix A. /// The order of the square matrix . - /// The B matrix. + /// On input the B matrix; on output the X matrix. /// This is equivalent to the GETRF and GETRS LAPACK routines. public virtual void LUSolve(int columnsOfB, float[] a, int order, float[] b) { @@ -1046,8 +1046,10 @@ namespace MathNet.Numerics.Algorithms.LinearAlgebra } var ipiv = new int[order]; - LUFactor(a, order, ipiv); - LUSolveFactored(columnsOfB, a, order, ipiv, b); + var clone = new float[a.Length]; + Buffer.BlockCopy(a, 0, clone, 0, a.Length * Constants.SizeOfFloat); + LUFactor(clone, order, ipiv); + LUSolveFactored(columnsOfB, clone, order, ipiv, b); } /// @@ -1057,7 +1059,7 @@ namespace MathNet.Numerics.Algorithms.LinearAlgebra /// The factored A matrix. /// The order of the square matrix . /// The pivot indices of . - /// The B matrix. + /// On input the B matrix; on output the X matrix. /// This is equivalent to the GETRS LAPACK routine. public virtual void LUSolveFactored(int columnsOfB, float[] a, int order, int[] ipiv, float[] b) { @@ -1146,103 +1148,6 @@ namespace MathNet.Numerics.Algorithms.LinearAlgebra } } - /// - /// Solves A*X=B for X using LU factorization. - /// - /// How to transpose the matrix. - /// The number of columns of B. - /// The square matrix A. - /// The order of the square matrix . - /// The B matrix. - /// This is equivalent to the GETRF and GETRS LAPACK routines. - public virtual void LUSolve(Transpose transposeA, int columnsOfB, float[] a, int order, float[] b) - { - if (a == null) - { - throw new ArgumentNullException("a"); - } - - if (b == null) - { - throw new ArgumentNullException("b"); - } - - if (a.Length != order * order) - { - throw new ArgumentException(Resources.ArgumentArraysSameLength, "a"); - } - - if (b.Length != order * columnsOfB) - { - throw new ArgumentException(Resources.ArgumentArraysSameLength, "b"); - } - - var ipiv = new int[order]; - LUFactor(a, order, ipiv); - LUSolveFactored(transposeA, columnsOfB, a, order, ipiv, b); - } - - /// - /// Solves A*X=B for X using a previously factored A matrix. - /// - /// How to transpose the matrix. - /// The number of columns of B. - /// The factored A matrix. - /// The order of the square matrix . - /// The pivot indices of . - /// The B matrix. - /// This is equivalent to the GETRS LAPACK routine. - public virtual void LUSolveFactored(Transpose transposeA, int columnsOfB, float[] a, int order, int[] ipiv, float[] b) - { - if (a == null) - { - throw new ArgumentNullException("a"); - } - - if (ipiv == null) - { - throw new ArgumentNullException("ipiv"); - } - - if (b == null) - { - throw new ArgumentNullException("b"); - } - - if (a.Length != order * order) - { - throw new ArgumentException(Resources.ArgumentArraysSameLength, "a"); - } - - if (ipiv.Length != order) - { - throw new ArgumentException(Resources.ArgumentArraysSameLength, "ipiv"); - } - - if (b.Length != order * columnsOfB) - { - throw new ArgumentException(Resources.ArgumentArraysSameLength, "b"); - } - - if ((transposeA == Transpose.Transpose) || (transposeA == Transpose.ConjugateTranspose)) - { - var aT = new float[a.Length]; - for (var i = 0; i < order; i++) - { - for (var j = 0; j < order; j++) - { - aT[(j * order) + i] = a[(i * order) + j]; - } - } - - LUSolveFactored(columnsOfB, aT, order, ipiv, b); - } - else - { - LUSolveFactored(columnsOfB, a, order, ipiv, b); - } - } - /// /// Computes the Cholesky factorization of A. /// diff --git a/src/Numerics/Algorithms/LinearAlgebra/native.common.include b/src/Numerics/Algorithms/LinearAlgebra/native.common.include index c12d8513..e59856ab 100644 --- a/src/Numerics/Algorithms/LinearAlgebra/native.common.include +++ b/src/Numerics/Algorithms/LinearAlgebra/native.common.include @@ -8,6 +8,7 @@ /// /// The requested of the matrix. /// + [SecuritySafeCritical] public override float MatrixNorm(Norm norm, int rows, int columns, float[] matrix) { if (matrix == null) @@ -46,6 +47,7 @@ /// /// The requested of the matrix. /// + [SecuritySafeCritical] public override float MatrixNorm(Norm norm, int rows, int columns, float[] matrix, float[] work) { if (matrix == null) @@ -86,6 +88,7 @@ /// /// The requested of the matrix. /// + [SecuritySafeCritical] public override double MatrixNorm(Norm norm, int rows, int columns, double[] matrix) { if (matrix == null) @@ -124,6 +127,7 @@ /// /// The requested of the matrix. /// + [SecuritySafeCritical] public override double MatrixNorm(Norm norm, int rows, int columns, double[] matrix, double[] work) { if (matrix == null) @@ -164,6 +168,7 @@ /// /// The requested of the matrix. /// + [SecuritySafeCritical] public override Complex32 MatrixNorm(Norm norm, int rows, int columns, Complex32[] matrix) { if (matrix == null) @@ -202,6 +207,7 @@ /// /// The requested of the matrix. /// + [SecuritySafeCritical] public override Complex32 MatrixNorm(Norm norm, int rows, int columns, Complex32[] matrix, float[] work) { if (matrix == null) @@ -242,6 +248,7 @@ /// /// The requested of the matrix. /// + [SecuritySafeCritical] public override Complex MatrixNorm(Norm norm, int rows, int columns, Complex[] matrix) { if (matrix == null) @@ -280,6 +287,7 @@ /// /// The requested of the matrix. /// + [SecuritySafeCritical] public override Complex MatrixNorm(Norm norm, int rows, int columns, Complex[] matrix, double[] work) { if (matrix == null) diff --git a/src/Numerics/Algorithms/LinearAlgebra/native.generic.include b/src/Numerics/Algorithms/LinearAlgebra/native.generic.include index 9a9d1f8e..b3f61c7c 100644 --- a/src/Numerics/Algorithms/LinearAlgebra/native.generic.include +++ b/src/Numerics/Algorithms/LinearAlgebra/native.generic.include @@ -6,6 +6,7 @@ /// The vector to add to . /// The result of the addition. /// This is similar to the AXPY BLAS routine. + [SecuritySafeCritical] public override void AddVectorToScaledVector(<#=dataType#>[] y, <#=dataType#> alpha, <#=dataType#>[] x, <#=dataType#>[] result) { if (y == null) @@ -43,6 +44,7 @@ /// The values to scale. /// This result of the scaling. /// This is similar to the SCAL BLAS routine. + [SecuritySafeCritical] public override void ScaleArray(<#=dataType#> alpha, <#=dataType#>[] x, <#=dataType#>[] result) { if (x == null) @@ -70,6 +72,7 @@ /// The vector y. /// The dot product of x and y. /// This is equivalent to the DOT BLAS routine. + [SecuritySafeCritical] public override <#=dataType#> DotProduct(<#=dataType#>[] x, <#=dataType#>[] y) { if (y == null) @@ -121,6 +124,7 @@ /// The number of columns in the matrix. /// The value to scale the matrix. /// The c matrix. + [SecuritySafeCritical] public override void MatrixMultiplyWithUpdate(Transpose transposeA, Transpose transposeB, <#=dataType#> alpha, <#=dataType#>[] a, int rowsA, int columnsA, <#=dataType#>[] b, int rowsB, int columnsB, <#=dataType#> beta, <#=dataType#>[] c) { if (a == null) @@ -164,6 +168,7 @@ /// The order of the square matrix . /// On exit, it contains the pivot indices. The size of the array must be . /// This is equivalent to the GETRF LAPACK routine. + [SecuritySafeCritical] public override void LUFactor(<#=dataType#>[] data, int order, int[] ipiv) { if (data == null) @@ -195,6 +200,7 @@ /// The N by N matrix to invert. Contains the inverse On exit. /// The order of the square matrix . /// This is equivalent to the GETRF and GETRI LAPACK routines. + [SecuritySafeCritical] public override void LUInverse(<#=dataType#>[] a, int order) { if (a == null) @@ -218,6 +224,7 @@ /// The order of the square matrix . /// The pivot indices of . /// This is equivalent to the GETRI LAPACK routine. + [SecuritySafeCritical] public override void LUInverseFactored(<#=dataType#>[] a, int order, int[] ipiv) { if (a == null) @@ -253,6 +260,7 @@ /// but should be N*blocksize. The blocksize is machine dependent. On exit, work[0] contains the optimal /// work size value. /// This is equivalent to the GETRF and GETRI LAPACK routines. + [SecuritySafeCritical] public override void LUInverse(<#=dataType#>[] a, int order, <#=dataType#>[] work) { if (a == null) @@ -288,6 +296,7 @@ /// but should be N*blocksize. The blocksize is machine dependent. On exit, work[0] contains the optimal /// work size value. /// This is equivalent to the GETRI LAPACK routine. + [SecuritySafeCritical] public override void LUInverseFactored(<#=dataType#>[] a, int order, int[] ipiv, <#=dataType#>[] work) { if (a == null) @@ -331,9 +340,25 @@ /// The order of the square matrix . /// The B matrix. /// This is equivalent to the GETRF and GETRS LAPACK routines. + [SecuritySafeCritical] public override void LUSolve(int columnsOfB, <#=dataType#>[] a, int order, <#=dataType#>[] b) { - throw new NotImplementedException(); + if (a == null) + { + throw new ArgumentNullException("a"); + } + + if (a.Length != order * order) + { + throw new ArgumentException(Resources.ArgumentArraysSameLength, "a"); + } + + if (b.Length != columnsOfB * order) + { + throw new ArgumentException(Resources.ArgumentArraysSameLength, "b"); + } + + SafeNativeMethods.<#=prefix#>_lu_solve(order, columnsOfB, a, b); } /// @@ -345,38 +370,35 @@ /// The pivot indices of . /// The B matrix. /// This is equivalent to the GETRS LAPACK routine. + [SecuritySafeCritical] public override void LUSolveFactored(int columnsOfB, <#=dataType#>[] a, int order, int[] ipiv, <#=dataType#>[] b) { - throw new NotImplementedException(); - } + if (a == null) + { + throw new ArgumentNullException("a"); + } - /// - /// Solves A*X=B for X using LU factorization. - /// - /// How to transpose the matrix. - /// The number of columns of B. - /// The square matrix A. - /// The order of the square matrix . - /// The B matrix. - /// This is equivalent to the GETRF and GETRS LAPACK routines. - public override void LUSolve(Transpose transposeA, int columnsOfB, <#=dataType#>[] a, int order, <#=dataType#>[] b) - { - throw new NotImplementedException(); - } + if (ipiv == null) + { + throw new ArgumentNullException("ipiv"); + } - /// - /// Solves A*X=B for X using a previously factored A matrix. - /// - /// How to transpose the matrix. - /// The number of columns of B. - /// The factored A matrix. - /// The order of the square matrix . - /// The pivot indices of . - /// The B matrix. - /// This is equivalent to the GETRS LAPACK routine. - public override void LUSolveFactored(Transpose transposeA, int columnsOfB, <#=dataType#>[] a, int order, int[] ipiv, <#=dataType#>[] b) - { - throw new NotImplementedException(); + if (a.Length != order * order) + { + throw new ArgumentException(Resources.ArgumentArraysSameLength, "a"); + } + + if (ipiv.Length != order) + { + throw new ArgumentException(Resources.ArgumentArraysSameLength, "ipiv"); + } + + if (b.Length != columnsOfB * order) + { + throw new ArgumentException(Resources.ArgumentArraysSameLength, "b"); + } + + SafeNativeMethods.<#=prefix#>_lu_solve_factored(order, columnsOfB, a, ipiv, b); } /// @@ -386,6 +408,7 @@ /// the Cholesky factorization. /// The number of rows or columns in the matrix. /// This is equivalent to the POTRF LAPACK routine. + [SecuritySafeCritical] public override void CholeskyFactor(<#=dataType#>[] a, int order) { if (a == null) @@ -416,6 +439,7 @@ /// The number of columns in the B matrix. /// This is equivalent to the POTRF add POTRS LAPACK routines. /// + [SecuritySafeCritical] public override void CholeskySolve(<#=dataType#>[] a, int orderA, <#=dataType#>[] b, int rowsB, int columnsB) { throw new NotImplementedException(); @@ -430,6 +454,7 @@ /// The number of rows in the B matrix. /// The number of columns in the B matrix. /// This is equivalent to the POTRS LAPACK routine. + [SecuritySafeCritical] public override void CholeskySolveFactored(<#=dataType#>[] a, int orderA, <#=dataType#>[] b, int rowsB, int columnsB) { throw new NotImplementedException(); @@ -445,6 +470,7 @@ /// On exit, A M by M matrix that holds the Q matrix of the /// QR factorization. /// This is similar to the GEQRF and ORGQR LAPACK routines. + [SecuritySafeCritical] public override void QRFactor(<#=dataType#>[] r, int rowsR, int columnsR, <#=dataType#>[] q) { throw new NotImplementedException(); @@ -463,6 +489,7 @@ /// but should be N*blocksize. The blocksize is machine dependent. On exit, work[0] contains the optimal /// work size value. /// This is similar to the GEQRF and ORGQR LAPACK routines. + [SecuritySafeCritical] public override void QRFactor(<#=dataType#>[] r, int rowsR, int columnsR, <#=dataType#>[] q, <#=dataType#>[] work) { throw new NotImplementedException(); @@ -480,6 +507,7 @@ /// The B matrix. /// The number of columns of B. /// On exit, the solution matrix. + [SecuritySafeCritical] public override void QRSolve(<#=dataType#>[] r, int rowsR, int columnsR, <#=dataType#>[] q, <#=dataType#>[] b, int columnsB, <#=dataType#>[] x) { throw new NotImplementedException(); @@ -500,6 +528,7 @@ /// The work array. The array must have a length of at least N, /// but should be N*blocksize. The blocksize is machine dependent. On exit, work[0] contains the optimal /// work size value. + [SecuritySafeCritical] public override void QRSolve(<#=dataType#>[] r, int rowsR, int columnsR, <#=dataType#>[] q, <#=dataType#>[] b, int columnsB, <#=dataType#>[] x, <#=dataType#>[] work) { throw new NotImplementedException(); @@ -515,6 +544,7 @@ /// The B matrix. /// The number of columns of B. /// On exit, the solution matrix. + [SecuritySafeCritical] public override void QRSolveFactored(<#=dataType#>[] q, <#=dataType#>[] r, int rowsR, int columnsR, <#=dataType#>[] b, int columnsB, <#=dataType#>[] x) { throw new NotImplementedException(); @@ -533,6 +563,7 @@ /// If is true, on exit VT contains the transposed /// right singular vectors. /// This is equivalent to the GESVD LAPACK routine. + [SecuritySafeCritical] public override void SingularValueDecomposition(bool computeVectors, <#=dataType#>[] a, int rowsA, int columnsA, <#=dataType#>[] s, <#=dataType#>[] u, <#=dataType#>[] vt) { throw new NotImplementedException(); @@ -554,6 +585,7 @@ /// Max(3*Min(M, N) + Max(M, N), 5*Min(M,N)). For complex matrices, 2*Min(M, N) + Max(M, N). /// On exit, work[0] contains the optimal work size value. /// This is equivalent to the GESVD LAPACK routine. + [SecuritySafeCritical] public override void SingularValueDecomposition(bool computeVectors, <#=dataType#>[] a, int rowsA, int columnsA, <#=dataType#>[] s, <#=dataType#>[] u, <#=dataType#>[] vt, <#=dataType#>[] work) { throw new NotImplementedException(); @@ -571,6 +603,7 @@ /// The B matrix. /// The number of columns of B. /// On exit, the solution matrix. + [SecuritySafeCritical] public override void SvdSolve(<#=dataType#>[] a, int rowsA, int columnsA, <#=dataType#>[] s, <#=dataType#>[] u, <#=dataType#>[] vt, <#=dataType#>[] b, int columnsB, <#=dataType#>[] x) { throw new NotImplementedException(); @@ -591,6 +624,7 @@ /// The work array. For real matrices, the work array should be at least /// Max(3*Min(M, N) + Max(M, N), 5*Min(M,N)). For complex matrices, 2*Min(M, N) + Max(M, N). /// On exit, work[0] contains the optimal work size value. + [SecuritySafeCritical] public override void SvdSolve(<#=dataType#>[] a, int rowsA, int columnsA, <#=dataType#>[] s, <#=dataType#>[] u, <#=dataType#>[] vt, <#=dataType#>[] b, int columnsB, <#=dataType#>[] x, <#=dataType#>[] work) { throw new NotImplementedException(); @@ -607,6 +641,7 @@ /// The B matrix. /// The number of columns of B. /// On exit, the solution matrix. + [SecuritySafeCritical] public override void SvdSolveFactored(int rowsA, int columnsA, <#=dataType#>[] s, <#=dataType#>[] u, <#=dataType#>[] vt, <#=dataType#>[] b, int columnsB, <#=dataType#>[] x) { throw new NotImplementedException(); diff --git a/src/Numerics/Algorithms/LinearAlgebra/native.header.include b/src/Numerics/Algorithms/LinearAlgebra/native.header.include index 3827a9aa..e14dfec2 100644 --- a/src/Numerics/Algorithms/LinearAlgebra/native.header.include +++ b/src/Numerics/Algorithms/LinearAlgebra/native.header.include @@ -36,6 +36,7 @@ namespace MathNet.Numerics.Algorithms.LinearAlgebra.<#=library#> { using System; using System.Numerics; + using System.Security; using Properties; /// diff --git a/src/Numerics/Algorithms/LinearAlgebra/safe.native.common.include b/src/Numerics/Algorithms/LinearAlgebra/safe.native.common.include index 230e00b8..4f72d5c6 100644 --- a/src/Numerics/Algorithms/LinearAlgebra/safe.native.common.include +++ b/src/Numerics/Algorithms/LinearAlgebra/safe.native.common.include @@ -40,6 +40,7 @@ namespace MathNet.Numerics.Algorithms.LinearAlgebra.<#= namespaceSuffix #> /// P/Invoke methods to the native math libraries. /// [SuppressUnmanagedCodeSecurity] + [SecurityCritical] internal static class SafeNativeMethods { /// @@ -161,27 +162,27 @@ namespace MathNet.Numerics.Algorithms.LinearAlgebra.<#= namespaceSuffix #> [DllImport(DllName, ExactSpelling = true, SetLastError = false, CallingConvention = CallingConvention.Cdecl)] internal static extern void z_lu_inverse_factored(int n, [In, Out] Complex[] a, [In, Out] int[] ipiv, [In, Out] Complex[] work, int lwork); - /* [DllImport(DllName, ExactSpelling = true, SetLastError = false, CallingConvention = CallingConvention.Cdecl)] - internal static extern void s_lu_solve_factored(int n, int nrhs, float a[], int ipiv[], float b[]) + [DllImport(DllName, ExactSpelling = true, SetLastError = false, CallingConvention = CallingConvention.Cdecl)] + internal static extern void s_lu_solve_factored(int n, int nrhs, float[] a, [In, Out]int[] ipiv, [In, Out] float[] b); [DllImport(DllName, ExactSpelling = true, SetLastError = false, CallingConvention = CallingConvention.Cdecl)] - internal static extern void d_lu_solve_factored(int n, int nrhs, double a[], int ipiv[], double b[]) + internal static extern void d_lu_solve_factored(int n, int nrhs, double[] a, [In, Out] int[] ipiv, [In, Out] double[] b); [DllImport(DllName, ExactSpelling = true, SetLastError = false, CallingConvention = CallingConvention.Cdecl)] - internal static extern void c_lu_solve_factored(int n, int nrhs, MKL_Complex8 a[], int ipiv[], MKL_Complex8 b[]) + internal static extern void c_lu_solve_factored(int n, int nrhs, Complex32[] a, [In, Out] int[] ipiv, [In, Out] Complex32[] b); [DllImport(DllName, ExactSpelling = true, SetLastError = false, CallingConvention = CallingConvention.Cdecl)] - internal static extern void z_lu_solve_factored(int n, int nrhs, MKL_Complex16 a[], int ipiv[], MKL_Complex16 b[]) + internal static extern void z_lu_solve_factored(int n, int nrhs, Complex[] a, [In, Out]int[] ipiv, [In, Out] Complex[] b); [DllImport(DllName, ExactSpelling = true, SetLastError = false, CallingConvention = CallingConvention.Cdecl)] - internal static extern void s_lu_solve(int n, int nrhs, float a[], float b[]) + internal static extern void s_lu_solve(int n, int nrhs, float[] a, [In, Out] float[] b); [DllImport(DllName, ExactSpelling = true, SetLastError = false, CallingConvention = CallingConvention.Cdecl)] - internal static extern void d_lu_solve(int n, int nrhs, double a[], double b[]) + internal static extern void d_lu_solve(int n, int nrhs, double[] a, [In, Out] double[] b); [DllImport(DllName, ExactSpelling = true, SetLastError = false, CallingConvention = CallingConvention.Cdecl)] - internal static extern void c_lu_solve(int n, int nrhs, MKL_Complex8 a[], MKL_Complex8 b[]) + internal static extern void c_lu_solve(int n, int nrhs, Complex32[] a, [In, Out] Complex32[] b); [DllImport(DllName, ExactSpelling = true, SetLastError = false, CallingConvention = CallingConvention.Cdecl)] - internal static extern void z_lu_solve(int n, int nrhs, MKL_Complex16 a[], MKL_Complex16 b[])*/ + internal static extern void z_lu_solve(int n, int nrhs, Complex[] a, [In, Out] Complex[] b); #endregion LAPACK diff --git a/src/Numerics/LinearAlgebra/Complex/Factorization/DenseLU.cs b/src/Numerics/LinearAlgebra/Complex/Factorization/DenseLU.cs index 0233dd39..fad99084 100644 --- a/src/Numerics/LinearAlgebra/Complex/Factorization/DenseLU.cs +++ b/src/Numerics/LinearAlgebra/Complex/Factorization/DenseLU.cs @@ -32,6 +32,7 @@ namespace MathNet.Numerics.LinearAlgebra.Complex.Factorization { using System; using System.Numerics; + using Algorithms.LinearAlgebra; using Generic; using Properties; using Threading; @@ -125,7 +126,7 @@ namespace MathNet.Numerics.LinearAlgebra.Complex.Factorization // LU solve by overwriting result. var dfactors = (DenseMatrix)Factors; - Control.LinearAlgebraProvider.LUSolveFactored(input.ColumnCount, dfactors.Data, dfactors.RowCount, Pivots, dresult.Data); + Control.LinearAlgebraProvider.LUSolveFactored(input.ColumnCount, dfactors.Data, dfactors.RowCount, Pivots, dresult.Data); } /// diff --git a/src/Numerics/LinearAlgebra/Complex32/Factorization/DenseLU.cs b/src/Numerics/LinearAlgebra/Complex32/Factorization/DenseLU.cs index 61ad70ed..ed68310f 100644 --- a/src/Numerics/LinearAlgebra/Complex32/Factorization/DenseLU.cs +++ b/src/Numerics/LinearAlgebra/Complex32/Factorization/DenseLU.cs @@ -31,6 +31,7 @@ namespace MathNet.Numerics.LinearAlgebra.Complex32.Factorization { using System; + using Algorithms.LinearAlgebra; using Generic; using Numerics; using Properties; @@ -125,7 +126,7 @@ namespace MathNet.Numerics.LinearAlgebra.Complex32.Factorization // LU solve by overwriting result. var dfactors = (DenseMatrix)Factors; - Control.LinearAlgebraProvider.LUSolveFactored(input.ColumnCount, dfactors.Data, dfactors.RowCount, Pivots, dresult.Data); + Control.LinearAlgebraProvider.LUSolveFactored(input.ColumnCount, dfactors.Data, dfactors.RowCount, Pivots, dresult.Data); } /// diff --git a/src/Numerics/LinearAlgebra/Double/Factorization/DenseLU.cs b/src/Numerics/LinearAlgebra/Double/Factorization/DenseLU.cs index 8b33dc87..78628204 100644 --- a/src/Numerics/LinearAlgebra/Double/Factorization/DenseLU.cs +++ b/src/Numerics/LinearAlgebra/Double/Factorization/DenseLU.cs @@ -31,6 +31,7 @@ namespace MathNet.Numerics.LinearAlgebra.Double.Factorization { using System; + using Algorithms.LinearAlgebra; using Generic; using Properties; diff --git a/src/Numerics/LinearAlgebra/Single/Factorization/DenseLU.cs b/src/Numerics/LinearAlgebra/Single/Factorization/DenseLU.cs index 8d199837..8b27c976 100644 --- a/src/Numerics/LinearAlgebra/Single/Factorization/DenseLU.cs +++ b/src/Numerics/LinearAlgebra/Single/Factorization/DenseLU.cs @@ -31,6 +31,7 @@ namespace MathNet.Numerics.LinearAlgebra.Single.Factorization { using System; + using Algorithms.LinearAlgebra; using Generic; using Properties; @@ -123,7 +124,7 @@ namespace MathNet.Numerics.LinearAlgebra.Single.Factorization // LU solve by overwriting result. var dfactors = (DenseMatrix)Factors; - Control.LinearAlgebraProvider.LUSolveFactored(input.ColumnCount, dfactors.Data, dfactors.RowCount, Pivots, dresult.Data); + Control.LinearAlgebraProvider.LUSolveFactored(input.ColumnCount, dfactors.Data, dfactors.RowCount, Pivots, dresult.Data); } /// diff --git a/src/UnitTests/LinearAlgebraProviderTests/Double/LinearAlgebraProviderTests.cs b/src/UnitTests/LinearAlgebraProviderTests/Double/LinearAlgebraProviderTests.cs index 39f48c5c..bc70778b 100644 --- a/src/UnitTests/LinearAlgebraProviderTests/Double/LinearAlgebraProviderTests.cs +++ b/src/UnitTests/LinearAlgebraProviderTests/Double/LinearAlgebraProviderTests.cs @@ -30,6 +30,9 @@ namespace MathNet.Numerics.UnitTests.LinearAlgebraProviderTests.Double using System.Collections.Generic; using Algorithms.LinearAlgebra; using LinearAlgebra.Double; + + using MathNet.Numerics.LinearAlgebra.Generic; + using NUnit.Framework; /// @@ -548,5 +551,70 @@ namespace MathNet.Numerics.UnitTests.LinearAlgebraProviderTests.Double AssertHelpers.AlmostEqual(a[8], -0.113636363636364, 14); } + /// + /// Can solve Ax=b using LU factorization. + /// + [Test] + public void CanSolveUsingLU() + { + var matrix = _matrices["Square3x3"]; + var a = new double[matrix.RowCount * matrix.RowCount]; + Array.Copy(matrix.Data, a, a.Length); + + var b = new[] { 1.0, 2.0, 3.0, 4.0, 5.0, 6.0 }; + Provider.LUSolve(2, a, matrix.RowCount, b); + + AssertHelpers.AlmostEqual(b[0], -1.477272727272726, 14); + AssertHelpers.AlmostEqual(b[1], -4.318181818181815, 14); + AssertHelpers.AlmostEqual(b[2], 3.068181818181816, 14); + AssertHelpers.AlmostEqual(b[3], -4.204545454545451, 14); + AssertHelpers.AlmostEqual(b[4], -12.499999999999989, 14); + AssertHelpers.AlmostEqual(b[5], 8.522727272727266, 14); + + NotModified(matrix.RowCount, matrix.ColumnCount, a, matrix); + } + + /// + /// Can solve Ax=b using LU factorization using a factored matrix. + /// + [Test] + public void CanSolveUsingLUOnFactoredMatrix() + { + var matrix = _matrices["Square3x3"]; + var a = new double[matrix.RowCount * matrix.RowCount]; + Array.Copy(matrix.Data, a, a.Length); + + var ipiv = new int[matrix.RowCount]; + Provider.LUFactor(a, matrix.RowCount, ipiv); + + var b = new[] { 1.0, 2.0, 3.0, 4.0, 5.0, 6.0 }; + Provider.LUSolveFactored(2, a, matrix.RowCount, ipiv, b); + + AssertHelpers.AlmostEqual(b[0], -1.477272727272726, 14); + AssertHelpers.AlmostEqual(b[1], -4.318181818181815, 14); + AssertHelpers.AlmostEqual(b[2], 3.068181818181816, 14); + AssertHelpers.AlmostEqual(b[3], -4.204545454545451, 14); + AssertHelpers.AlmostEqual(b[4], -12.499999999999989, 14); + AssertHelpers.AlmostEqual(b[5], 8.522727272727266, 14); + } + + /// + /// Checks to see if a matrix and array contain the same values. + /// + /// number of rows. + /// number of columns. + /// array to check. + /// matrix to check against. + private static void NotModified(int rows, int columns, IList array, Matrix matrix) + { + var index = 0; + for (var col = 0; col < columns; col++) + { + for (var row = 0; row < rows; row++) + { + Assert.AreEqual(array[index++], matrix[row, col]); + } + } + } } }