DenseMat.cpp

Go to the documentation of this file.

 
#include "DenseMat.hpp"
 
void MinEigenSY(const int& N, float* A, float* w, float* work, const int& lwork)
{
    int info;
    int iwork[1] = { 0 };
    int liwork = 1;
    char uplo{ 'U' };
    char Nonly{ 'N' };
 
    ssyevd_(&Nonly, &uplo, &N, A, &N, w, work, &lwork, iwork, &liwork, &info);
    if (info > 0) {
        cout << "failed to compute eigenvalues!" << endl;
    }
}
 
 
//void MinEigenS(const int& N, float* a, float* w)
//{
//    MKL_INT info = LAPACKE_ssyevd(LAPACK_COL_MAJOR, 'N', 'U', N, a, N, w);
//    if (info > 0) {
//        cout << "failed to compute eigenvalues!" << endl;
//    }
//}
 
void Dcopy(const int& N, double* dst, const double* src)
{
    const int incx = 1, incy = 1;
    dcopy_(&N, src, &incx, dst, &incy);
}
 
double Ddot(int n, double* a, double* b)
{
    const int inca = 1, incb = 1;
    return ddot_(&n, a, &inca, b, &incb);
}
 
// WARNING: absolute sum!
double Dnorm1(const int& N, double* x)
{
    const int incx = 1;
    return dasum_(&N, x, &incx);
}
 
double Dnorm2(const int& N, double* x)
{
    const int incx = 1;
    return dnrm2_(&N, x, &incx);
}
 
void Dscalar(const int& n, const double& alpha, double* x)
{
    // x = a x
    const int incx = 1;
    dscal_(&n, &alpha, x, &incx);
}
 
void Daxpy(const int& n, const double& alpha, const double* x, double* y)
{
    // y= ax +y
    const int incx = 1, incy = 1;
    daxpy_(&n, &alpha, x, &incx, y, &incy);
}
 
void DaABpbC(const int& m, const int& n, const int& k, const double& alpha,
             const double* A, const double* B, const double& beta, double* C)
{
    /*  C' = alpha B'A' + beta C'
     *  A: m x k
     *  B: k x n
     *  C: m x n
     *  all column majored matrices, no tranpose
     *  A' in col-order in Fortran = A in row-order in C/Cpp
     */
 
    const char transa = 'N', transb = 'N';
    dgemm_(&transa, &transb, &n, &m, &k, &alpha, B, &n, A, &k, &beta, C, &n);
}
 
 
void myDABpC(const int& m, const int& n, const int& k, const double* A, const double* B, double* C)
{
    // C = AB + C
    // A: m*n  B:n*k  C:m*k
    // all matrix are row majored matrices
 
    for (int i = 0; i < m; i++) {
        for (int j = 0; j < k; j++) {
            for (int m = 0; m < n; m++) {
                C[i * k + j] += A[i * n + m] * B[m * k + j];
            }
        }
    }
}
 
void myDABpCp(const int& m, const int& n, const int& k, const double* A, const double* B, double* C, const int* flag, const int N)
{
    for (int i = 0; i < m; i++) {
        for (int j = 0; j < k; j++) {
            for (int p = 0; p < 3; p++) {
                if (flag[p] != 0)
                    C[i * k + j] += A[i * n + p] * B[p * k + j];
            }
            for (int p = 0; p < 2; p++) {
                if (flag[p] != 0) {
                    for (int m = 0; m < N; m++) {
                        C[i * k + j] += A[i * n + 3 + p * (N + 1) + m] * B[(3 + p * (N + 1) + m) * k + j];
                    }
                }
            }
        }
    }
}
 
void myDABpCp1(const int& m, const int& n, const int& k, const double* A, const double* B, double* C, const int* flag, const int N)
{
 
    for (int i = 0; i < m; i++) {
        for (int j = 0; j < k; j++) {
            int s = 0;
            for (int p = 0; p < 3; p++) {
                if (flag[p] != 0) {
                    C[i * k + j] += A[i * n + p] * B[s * k + j];
                    s++;
                }                   
            }
            for (int p = 0; p < 2; p++) {
                if (flag[p] != 0) {
                    for (int m = 0; m < N; m++) {
                        C[i * k + j] += A[i * n + 3 + p * (N + 1) + m] * B[s * k + j];
                        s++;
                    }
                }
            }
        }
    }
}
 
void myDABpCp2(const int& m, const int& n, const int& k, const double* A, const double* B, double* C, const int* flag, const int N)
{
    for (int i = 0; i < m; i++) {
        for (int j = 0; j < k; j++) {
            int s = 0;
            for (int p = 0; p < 3; p++) {
                if (flag[p] != 0) {
                    C[i * k + j] += A[i * n + s] * B[p * k + j];
                    s++;
                }                   
            }
            for (int p = 0; p < 2; p++) {
                if (flag[p] != 0) {
                    for (int m = 0; m < N; m++) {
                        C[i * k + j] += A[i * n + s] * B[(3 + p * (N + 1) + m) * k + j];
                        s++;
                    }
                }
            }
        }
    }
}
 
void DaAxpby(const int& m, const int& n, const double& a, const double* A,
             const double* x, const double& b, double* y)
{
    /*  y= aAx+by
     */
    for (int i = 0; i < m; i++) {
        y[i] = b * y[i];
        for (int j = 0; j < n; j++) {
            y[i] += a * A[i * n + j] * x[j];
        }
    }
}
 
void LUSolve(const int& nrhs, const int& N, double* A, double* b, int* pivot)
{
    int info;
 
    dgesv_(&N, &nrhs, A, &N, pivot, b, &N, &info);
 
    if (info < 0) {
        cout << "Wrong Input !" << endl;
    } else if (info > 0) {
        cout << "Singular Matrix !" << endl;
    }
}
 
void SYSSolve(const int& nrhs, const char* uplo, const int& N, double* A, double* b, int* pivot, double* work, const int& lwork)
{
    int    info;
 
    dsysv_(uplo, &N, &nrhs, A, &N, pivot, b, &N, work, &lwork, &info);
    if (info < 0) {
        cout << "Wrong Input !" << endl;
    } else if (info > 0) {
        cout << "Singular Matrix !" << endl;
    }
}
 
 
/*----------------------------------------------------------------------------*/
/*  Brief Change History of This File                                         */
/*----------------------------------------------------------------------------*/
/*  Author              Date             Actions                              */
/*----------------------------------------------------------------------------*/
/*  Shizhe Li           Oct/21/2021      Create file                          */
/*----------------------------------------------------------------------------*/