Invers matrice

Hej

Prøv at kigge om der ikke findes en funktion i LAPACK (http://www.netlib.org/lapack/)

LAPACK kan goere det paa den hele rigtige maade. Men er en "stor pakke".

Matrix inversion kan laves rimeligt simpelt hvis man ikke har store krav til numerisk stabilitet.

Her er noget C kode (som er omskrevet fra noget C++ kode som jeg lavede for mange aar siden):

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

#define IX(row,col,ncol) (ncol*row+col)

void print(double *a, int nrow, int ncol)
{
    int i, j;
    for(i = 0; i < nrow; i++)
    {
        for(j = 0; j < ncol; j++)
        {
            printf(" %9.2f", a[IX(i,j,ncol)]);
        }
        printf("\n");
    }
}

/* C = A*B */
void mul(double *a, double *b, double *c, int nrow, int ncol, int dim3)
{
    int i, j, k;
    double tmp;
    for(i = 0; i < nrow; i++)
    {
        for(j = 0; j < ncol; j++) {
            tmp = 0;
            for(k = 0; k < dim3; k++) tmp = tmp + a[IX(i,k,dim3)] * b[IX(k,j,ncol)];
            c[IX(i,j,ncol)] = tmp;
        }
    }
}

void submulrow(double *x, int ncol, int r1, int r2, double scale)
{
    int j;
    for(j = 0; j < ncol; j++) x[IX(r1,j,ncol)] =  x[IX(r1,j,ncol)] - scale*x[IX(r2,j,ncol)];
    return;
}

/* B = INV(A) using Gauss */
void inv(double *a, double* b, int nrow, int ncol)
{
    int i, j, k;
    double *c;
    double tmp;
    c = (double *)malloc(2*nrow*ncol*sizeof(double));
    for(i = 0; i < nrow; i++)
    {
        for(j = 0; j < ncol; j++) c[IX(i,j,2*ncol)] = a[IX(i,j,ncol)];
        for(j = ncol; j < 2*ncol; j++) c[IX(i,j,2*ncol)] = (i+ncol==j) ? 1 : 0;
    }
    for(i = 0; i < nrow; i++)
    {
        for(k = (i+1); k < nrow; k++) if(c[IX(i,i,2*ncol)] == 0) submulrow(c, 2*ncol, i, (i+k)%nrow, 1.0);
        for(k = 0; k < ncol; k++)
        {
            if(i != k)
            {
                tmp = c[IX(k,i,2*ncol)] / c[IX(i,i,2*ncol)];
                submulrow(c, 2*ncol, k, i, tmp);
            }
            else
            {
                tmp = c[IX(k,k,2*ncol)];
                for(j = 0; j < 2*ncol; j++) c[IX(k,j,2*ncol)] = c[IX(k,j,2*ncol)] / tmp;
            }
        }
    }
    for(i = 0; i < nrow; i++)
    {
        for(j = 0; j < ncol; j++) b[IX(i,j,ncol)] = c[IX(i,j+ncol,2*ncol)];
    }
    free(c);
}

int main()
{
    double x[2*2] = { 1, 2, 3, 4 };
    double xx[2*2];
    double invx[2*2];
    double xinvx[2*2];
    print(x, 2, 2);
    mul(x, x, xx, 2, 2, 2);
    print(xx, 2, 2);
    inv(x, invx, 2, 2);
    print(invx, 2, 2);
    mul(x, invx, xinvx, 2, 2, 2);
    print(xinvx, 2, 2);
    return EXIT_SUCCESS;
}


Den er ihvertfald rimelig simpel.

Jeg kigger på LAPACK og din kode i aften arne_v - tak så langt !

Eherm ... Det dudder ikke, men har fundet en løsning ... Jeg har været så smart at forsøge et invertere en 1x1 - altså - skalar ... Den simplestemåde at gøre deet på er vidst ved at tage en reciprokke værdi :)

Smid lige et svar begge to som tak for hjælpen ..

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

#define IX(row,col,ncol) (ncol*row+col)

void print(double *a, int nrow, int ncol)
{
    int i, j;
    for(i = 0; i < nrow; i++)
    {
        for(j = 0; j < ncol; j++)
        {
            printf(" %9.2f", a[IX(i,j,ncol)]);
        }
        printf("\n");
    }
}

/* C = A*B */
void mul(double *a, double *b, double *c, int nrow, int ncol, int dim3)
{
    int i, j, k;
    double tmp;
    for(i = 0; i < nrow; i++)
    {
        for(j = 0; j < ncol; j++) {
            tmp = 0;
            for(k = 0; k < dim3; k++) tmp = tmp + a[IX(i,k,dim3)] * b[IX(k,j,ncol)];
            c[IX(i,j,ncol)] = tmp;
        }
    }
}

void submulrow(double *x, int ncol, int r1, int r2, double scale)
{
    int j;
    for(j = 0; j < ncol; j++) x[IX(r1,j,ncol)] =  x[IX(r1,j,ncol)] - scale*x[IX(r2,j,ncol)];
    return;
}

/* B = INV(A) using Gauss */
void inv(double *a, double* b, int nrow, int ncol)
{
    int i, j, k;
    double *c;
    double tmp;
    c = (double *)malloc(2*nrow*ncol*sizeof(double));
    for(i = 0; i < nrow; i++)
    {
        for(j = 0; j < ncol; j++) c[IX(i,j,2*ncol)] = a[IX(i,j,ncol)];
        for(j = ncol; j < 2*ncol; j++) c[IX(i,j,2*ncol)] = (i+ncol==j) ? 1 : 0;
    }
    for(i = 0; i < nrow; i++)
    {
        for(k = (i+1); k < nrow; k++) if(c[IX(i,i,2*ncol)] == 0) submulrow(c, 2*ncol, i, (i+k)%nrow, 1.0);
        for(k = 0; k < ncol; k++)
        {
            if(i != k)
            {
                tmp = c[IX(k,i,2*ncol)] / c[IX(i,i,2*ncol)];
                submulrow(c, 2*ncol, k, i, tmp);
            }
            else
            {
                tmp = c[IX(k,k,2*ncol)];
                for(j = 0; j < 2*ncol; j++) c[IX(k,j,2*ncol)] = c[IX(k,j,2*ncol)] / tmp;
            }
        }
    }
    for(i = 0; i < nrow; i++)
    {
        for(j = 0; j < ncol; j++) b[IX(i,j,ncol)] = c[IX(i,j+ncol,2*ncol)];
    }
    free(c);
}

int main()
{
    double x[1*1] = { 2 };
    double invx[1*1];
    double xinvx[1*1];
    print(x, 1, 1);
    inv(x, invx, 1, 1);
    print(invx, 1, 1);
    mul(x, invx, xinvx, 1, 1, 1);
    print(xinvx, 1, 1);
    return EXIT_SUCCESS;
}

giver da ellers korrekt resultat hos mig.

Men anyway - et svar.