example_dpotrf.c File Reference

Example of Cholesky factorization. More...

#include <stdlib.h>
#include <stdio.h>
#include <string.h>
#include <math.h>
#include <plasma.h>
#include <cblas.h>
#include <lapacke.h>
#include <core_blas.h>

Include dependency graph for example_dpotrf.c:

Go to the source code of this file.

Functions
int	check_factorization (int, double *, double *, int, int)
int	main ()

Variables
int	IONE = 1
int	ISEED [4] = {0,0,0,1}

---

Detailed Description

Example of Cholesky factorization.

PLASMA testing routines PLASMA is a software package provided by Univ. of Tennessee, Univ. of California Berkeley and Univ. of Colorado Denver

Version:

2.4.5

Author:

Bilel Hadri

Date:

2010-11-15 d Tue Nov 22 14:35:53 2011

Definition in file example_dpotrf.c.

---

Function Documentation

int check_factorization	(	int	N,
		double *	A1,
		double *	A2,
		int	LDA,
		int	uplo
	)

Definition at line 80 of file example_dpotrf.c.

References cblas_dtrmm(), CblasColMajor, CblasLeft, CblasLower, CblasNonUnit, CblasRight, CblasTrans, CblasUpper, lapack_const, PlasmaInfNorm, and PlasmaUpper.

{
    double Anorm, Rnorm;
    double alpha;
    int info_factorization;
    int i,j;
    double eps;

    eps = LAPACKE_dlamch_work('e');

    double *Residual = (double *)malloc(N*N*sizeof(double));
    double *L1       = (double *)malloc(N*N*sizeof(double));
    double *L2       = (double *)malloc(N*N*sizeof(double));
    double *work              = (double *)malloc(N*sizeof(double));

    memset((void*)L1, 0, N*N*sizeof(double));
    memset((void*)L2, 0, N*N*sizeof(double));

    alpha= 1.0;

    LAPACKE_dlacpy_work(LAPACK_COL_MAJOR,' ', N, N, A1, LDA, Residual, N);

    /* Dealing with L'L or U'U  */
    if (uplo == PlasmaUpper){
        LAPACKE_dlacpy_work(LAPACK_COL_MAJOR,'u', N, N, A2, LDA, L1, N);
        LAPACKE_dlacpy_work(LAPACK_COL_MAJOR,'u', N, N, A2, LDA, L2, N);
        cblas_dtrmm(CblasColMajor, CblasLeft, CblasUpper, CblasTrans, CblasNonUnit, N, N, (alpha), L1, N, L2, N);
    }
    else{
        LAPACKE_dlacpy_work(LAPACK_COL_MAJOR,'l', N, N, A2, LDA, L1, N);
        LAPACKE_dlacpy_work(LAPACK_COL_MAJOR,'l', N, N, A2, LDA, L2, N);
        cblas_dtrmm(CblasColMajor, CblasRight, CblasLower, CblasTrans, CblasNonUnit, N, N, (alpha), L1, N, L2, N);
    }

    /* Compute the Residual || A -L'L|| */
    for (i = 0; i < N; i++)
        for (j = 0; j < N; j++)
           Residual[j*N+i] = L2[j*N+i] - Residual[j*N+i];

    Rnorm = LAPACKE_dlange_work(LAPACK_COL_MAJOR, lapack_const(PlasmaInfNorm), N, N, Residual, N, work);
    Anorm = LAPACKE_dlange_work(LAPACK_COL_MAJOR, lapack_const(PlasmaInfNorm), N, N, A1, LDA, work);

    printf("============\n");
    printf("Checking the Cholesky Factorization \n");
    printf("-- ||L'L-A||_oo/(||A||_oo.N.eps) = %e \n",Rnorm/(Anorm*N*eps));

    if ( isnan(Rnorm/(Anorm*N*eps)) || (Rnorm/(Anorm*N*eps) > 10.0) ){
        printf("-- Factorization is suspicious ! \n");
        info_factorization = 1;
    }
    else{
        printf("-- Factorization is CORRECT ! \n");
        info_factorization = 0;
    }

    free(Residual); free(L1); free(L2); free(work);

    return info_factorization;
}

Here is the call graph for this function:

int main

(

)

Definition at line 32 of file example_dpotrf.c.

References check_factorization(), PLASMA_dlacpy, PLASMA_dplgsy(), PLASMA_dpotrf(), PLASMA_Finalize(), PLASMA_Init(), PlasmaUpper, and PlasmaUpperLower.

{

    int cores = 2;
    int N     = 10 ;
    int LDA   = 10 ;
    int info_factorization;

    double *A1   = (double *)malloc(LDA*N*sizeof(double));
    double *A2   = (double *)malloc(LDA*N*sizeof(double));

    /* Check if unable to allocate memory */
    if ((!A1)||(!A2)){
        printf("Out of Memory \n ");
        return EXIT_SUCCESS;
    }

    /* Plasma Initialize */
    PLASMA_Init(cores);
    printf("-- PLASMA is initialized to run on %d cores. \n",cores);

    /* Initialize A1 and A2 for Symmetric Positive Matrix */
    PLASMA_dplgsy( (double)N, N, A1, LDA, 51 );
    PLASMA_dlacpy( PlasmaUpperLower, N, N, A1, LDA, A2, LDA );

    /* Plasma routines */
    PLASMA_dpotrf(PlasmaUpper, N, A2, LDA);

    /* Check the factorization */
    info_factorization = check_factorization( N, A1, A2, LDA, PlasmaUpper);

    if ( info_factorization != 0 )
       printf("-- Error in DPOTRF example ! \n");
    else
       printf("-- Run of DPOTRF example successful ! \n");

    free(A1); free(A2);

    PLASMA_Finalize();

    return EXIT_SUCCESS;
}

Here is the call graph for this function:

---

Variable Documentation

int IONE = 1

Definition at line 29 of file example_dpotrf.c.

int ISEED[4] = {0,0,0,1}

Definition at line 30 of file example_dpotrf.c.