lapack.c

Go to the documentation of this file.

/*
 * This work is licensed under a Creative Commons 
 * Attribution-Noncommercial-Share Alike 3.0 United States License.
 * 
 *    http://creativecommons.org/licenses/by-nc-sa/3.0/us/
 * 
 * You are free:
 * 
 *    to Share - to copy, distribute, display, and perform the work
 *    to Remix - to make derivative works
 * 
 * Under the following conditions:
 * 
 *    Attribution. You must attribute the work in the manner specified by the
 *    author or licensor (but not in any way that suggests that they endorse you
 *    or your use of the work).
 * 
 *    Noncommercial. You may not use this work for commercial purposes.
 * 
 *    Share Alike. If you alter, transform, or build upon this work, you may
 *    distribute the resulting work only under the same or similar license to
 *    this one.
 * 
 * For any reuse or distribution, you must make clear to others the license
 * terms of this work. The best way to do this is by including this header.
 * 
 * Any of the above conditions can be waived if you get permission from the
 * copyright holder.
 * 
 * Apart from the remix rights granted under this license, nothing in this
 * license impairs or restricts the author's moral rights.
 */


#include <jwsc/config.h>

#include <stdlib.h>
#include <assert.h>

#include "jwsc/math/lapack.h"

#ifdef JWSC_HAVE_DMALLOC
#include <dmalloc.h>
#endif


extern float  sgetrf_(int*,int*,float*,int*,int*,int*);
extern double dgetrf_(int*,int*,double*,int*,int*,int*);

extern float  sgetri_(int*,float*,int*,int*,float*,int*,int*);
extern double dgetri_(int*,double*,int*,int*,double*,int*,int*);

extern float sgesdd_(char*,int*,int*,float*,int*,float*,float*,int*,float*,int*,float*,int*,int*,int*);
extern double dgesdd_(char*,int*,int*,double*,int*,double*,double*,int*,double*,int*,double*,int*,int*,int*);

extern float ssyev_(char*,char*,int*,float*,int*,float*,float*,int*,int*);
extern double dsyev_(char*,char*,int*,double*,int*,double*,double*,int*,int*);


static void trans_d(double** A_out, const double* A, int M, int N)
{
    int m, n;

    assert(*A_out != A);

    if (*A_out == NULL)
    {
        assert(*A_out = malloc(M*N*sizeof(double)));
    }

    for (m = 0; m < M; m++)
    {
        for (n = 0; n < N; n++)
        {
            (*A_out)[ n*M + m ] = A[ m*N + n ];
        }
    }
}


static void trans_f(float** A_out, const float* A, int M, int N)
{
    int i, j;

    assert(*A_out != A);

    if (*A_out == NULL)
    {
        assert(*A_out = malloc(M*N*sizeof(float)));
    }

    for (i = 0; i < M; i++)
    {
        for (j = 0; j < N; j++)
        {
            (*A_out)[ j*M + i ] = A[ i*N + j ];
        }
    }
}


void lapack_sgetrf
(
    int     M, 
    int     N, 
    float*  A, 
    int     lda,
    int*    ipiv,
    int*    info
)
{
    float*  A_t = NULL;

    trans_f(&A_t, A, M, N);
    sgetrf_(&N, &M, A_t, &lda, ipiv, info);
    trans_f(&A, A_t, N, M);

    free(A_t);
}


void lapack_dgetrf
(
    int     M, 
    int     N, 
    double* A, 
    int     lda,
    int*    ipiv,
    int*    info
)
{
    double* A_t = NULL;

    trans_d(&A_t, A, M, N);
    dgetrf_(&M, &N, A_t, &lda, ipiv, info);
    trans_d(&A, A_t, N, M);

    free(A_t);
}


void lapack_sgetri
(
    int     N, 
    float*  A, 
    int     lda,
    int*    ipiv,
    float*  work, 
    int     lwork,
    int*    info
)
{
    float*  A_t = NULL;

    trans_f(&A_t, A, N, N);
    sgetri_(&N, A_t, &lda, ipiv, work, &lwork, info);
    trans_f(&A, A_t, N, N);

    free(A_t);
}


void lapack_dgetri
(
    int     N, 
    double* A, 
    int     lda,
    int*    ipiv,
    double* work, 
    int     lwork,
    int*    info
)
{
    double* A_t = NULL;

    trans_d(&A_t, A, N, N);
    dgetri_(&N, A_t, &lda, ipiv, work, &lwork, info);
    trans_d(&A, A_t, N, N);

    free(A_t);
}


void lapack_sgesdd
(
    char   JOBZ,
    int    M,
    int    N, 
    float* A, 
    int    LDA,
    float* S,
    float* U,
    int    LDU,
    float* VT,
    int    LDVT,
    float* WORK,
    int    LWORK,
    int*   IWORK,
    int*   INFO
)
{
    int    min_MN;
    float* A_t  = NULL;
    float* U_t  = NULL;
    float* VT_t = NULL;

    min_MN = (M <= N) ? M : N;

    if (JOBZ == 'A')
    {
        trans_f(&U_t, U, LDU, M);
        trans_f(&VT_t, VT, LDVT, N);
    }
    else if (JOBZ == 'S')
    {
        trans_f(&U_t, U, LDU, min_MN);
        trans_f(&VT_t, VT, LDVT, min_MN);
    }
    else if (JOBZ == 'O')
    {
        if (M < N) 
            trans_f(&U_t, U, LDU, M);
        else
            trans_f(&VT_t, VT, LDVT, N);
    }
    else if (JOBZ == 'N')
    {
        ;
    }
    else
    {
        *INFO = -1;
        return;
    }

    trans_f(&A_t, A, M, N);
    sgesdd_(&JOBZ, &M, &N, A_t, &LDA, S, U_t, &LDU, VT_t, &LDVT, WORK,
            &LWORK, IWORK, INFO);
    trans_f(&A, A_t, N, M);

    if (JOBZ == 'A')
    {
        trans_f(&U, U_t, M, LDU);
        trans_f(&VT, VT_t, N, LDVT);
    }
    else if (JOBZ == 'S')
    {
        trans_f(&U, U_t, min_MN, LDU);
        trans_f(&VT, VT_t, min_MN, LDVT);
    }
    else if (JOBZ == 'O')
    {
        if (M < N) 
            trans_f(&U, U_t, M, LDU);
        else
            trans_f(&VT, VT_t, N, LDVT);
    }

    free(A_t);
    free(U_t);
    free(VT_t);
}


void lapack_dgesdd
(
    char    JOBZ,
    int     M,
    int     N, 
    double* A, 
    int     LDA,
    double* S,
    double* U,
    int     LDU,
    double* VT,
    int     LDVT,
    double* WORK,
    int     LWORK,
    int*    IWORK,
    int*    INFO
)
{
    int     min_MN;
    double* A_t  = NULL;
    double* U_t  = NULL;
    double* VT_t = NULL;

    min_MN = (M <= N) ? M : N;

    if (JOBZ == 'A')
    {
        trans_d(&U_t, U, LDU, M);
        trans_d(&VT_t, VT, LDVT, N);
    }
    else if (JOBZ == 'S')
    {
        trans_d(&U_t, U, LDU, min_MN);
        trans_d(&VT_t, VT, LDVT, min_MN);
    }
    else if (JOBZ == 'O')
    {
        if (M < N) 
            trans_d(&U_t, U, LDU, M);
        else
            trans_d(&VT_t, VT, LDVT, N);
    }
    else if (JOBZ == 'N')
    {
        ;
    }
    else
    {
        *INFO = -1;
        return;
    }

    trans_d(&A_t, A, M, N);
    dgesdd_(&JOBZ, &M, &N, A_t, &LDA, S, U_t, &LDU, VT_t, &LDVT, WORK,
            &LWORK, IWORK, INFO);
    trans_d(&A, A_t, N, M);

    if (JOBZ == 'A')
    {
        trans_d(&U, U_t, M, LDU);
        trans_d(&VT, VT_t, N, LDVT);
    }
    else if (JOBZ == 'S')
    {
        trans_d(&U, U_t, min_MN, LDU);
        trans_d(&VT, VT_t, min_MN, LDVT);
    }
    else if (JOBZ == 'O')
    {
        if (M < N) 
            trans_d(&U, U_t, M, LDU);
        else
            trans_d(&VT, VT_t, N, LDVT);
    }

    free(A_t);
    free(U_t);
    free(VT_t);
}


void lapack_ssyev
(
    char    JOBZ,
    char    UPLO,
    int     N, 
    float*  A, 
    int     LDA,
    float*  W,
    float*  WORK,
    int     LWORK,
    int*    INFO
)
{
    float* A_t = NULL;

    if (JOBZ != 'N' && JOBZ != 'V')
    {
        *INFO = -1;
        return;
    }

    trans_f(&A_t, A, N, N);
    ssyev_(&JOBZ, &UPLO, &N, A_t, &LDA, W, WORK, &LWORK, INFO);
    trans_f(&A, A_t, N, N);

    free(A_t);
}


void lapack_dsyev
(
    char    JOBZ,
    char    UPLO,
    int     N, 
    double* A, 
    int     LDA,
    double* W,
    double* WORK,
    int     LWORK,
    int*    INFO
)
{
    double* A_t = NULL;

    if (JOBZ != 'N' && JOBZ != 'V')
    {
        *INFO = -1;
        return;
    }

    trans_d(&A_t, A, N, N);
    dsyev_(&JOBZ, &UPLO, &N, A_t, &LDA, W, WORK, &LWORK, INFO);
    trans_d(&A, A_t, N, N);

    free(A_t);
}