Geek asks and answers

Solution of a linear system with dpotrs (Cholesky factorization)

I am trying to solve a linear system of equations with a valve.

My code is as follows:

//ATTENTION: matrix in column-major
double A[3*3]={  2.0, -1.0,  0.0,
                0.0,  2.0, -1.0,
                0.0,  0.0,  2.0},
b[3]={1.0,2.0,3.0};

integer n=3,info,nrhs=1; char uplo='L';

dpotrf_("L", &n, A, &n, &info);
dpotrs_("L", &n, &nrhs, A, &n, b, &n, &info);

printf("Solution: %10.4f %10.4f %10.4f",b[0], b[1], b[2]);

As asked in this question , you must first expand the matrix. dpotrf is supposed to factorize, dpotrs solves the system using a factorized matrix.

However my result

2.5   4.0   3.5

clearly wrong, I checked it here with WolframAlpha

Where is my mistake?

+3
source share
1 answer

Curious what 2.5 4.0 3.5is the rhs solution 1 2 3... if the matrix

 2   -1  0
 -1  2   -1
 0   -1  2

dpotrfand are dpotrsmade for symmetric positive definite matrices ...

I would suggest using dgetrfand dgetrs. In your case:

#include <stdio.h>
#include <string.h>
#include <stdlib.h>
#include <lapacke.h>

int main()
{
    //ATTENTION: matrix in column-major
    double A[3*3]={  2.0, -1.0,  0.0,
            0,  2.0, -1.0,
            0.0,  0,  2.0},
            b[3]={1.0,2,3.0};

    int n=3,info,nrhs=1; char uplo='L';
    int ipvs[3];
    dgetrf_(&n, &n, A, &n,ipvs, &info);
    printf("info %d\n",info);
    dgetrs_("T", &n, &nrhs, A, &n,ipvs, b, &n, &info);
    printf("info %d\n",info);
    printf("Solution: %10.4f %10.4f %10.4f\n",b[0], b[1], b[2]);

    return 0;
}

1.3750 1.75 1.5. , "N" "T". 0.5 1.25 2.125

!

,

+2

More articles:


All Articles