luwerwer.txt

通过LU分解求矩阵的逆矩阵

// LUqiuni.cpp : Defines the entry point for the console application.
//

#include "stdafx.h"
#include<math.h>
#include<stdio.h>
#include<iomanip.h>
#include <stdlib.h>

#define n 10

int LU(double a[n][n])
{
    int i,j,k,p;
	double b[n][n],u[n][n],l[n][n];
	double s;
	for(k=0;k<n;k++) /* LU分解 */
	{
		for(j=k;j<n;j++)
        {
			s=0.0;
		    for(p=0;p<=k-1;p++)
				s=s+a[k][p]*a[p][j];
		    a[k][j]=a[k][j]-s;
		}
		if(a[k][k]==0) {printf ("A not LU.\n"); return 0;}
		for(i=k+1;i<n;i++)
	    {
			s=0;
		    for(p=0;p<=k-1;p++)
				s=s+a[i][p]*a[p][k];
		    a[i][k]=a[i][k]-s/a[k][k];
		}
	}
	for (i=0;i<n;i++) /*求矩阵U的逆 */
    {
		u[i][i]=1/a[i][i];
		for(k=i-1;k>=0;k--)
        {
			s=0;
            for(j=k+1;j<=i;j++)
				s=s+a[k][j]*u[j][i];
            u[k][i]=-s/a[k][k];
		}
	}
	for(i=0;i<n;i++) /*求矩阵L的逆 */
    {
		l[i][i]=1;
		for(k=i+1;k<n;k++)
		{
			for(j=i;j<=k-1;j++)
				l[k][i]=l[k][i]-a[k][j]*l[j][i];
		}
	}
	for(i=0;i<n;i++) /* 计算A的逆 */
	{
		for(j=0;j<n;j++)
		{
			for(k=0;k<n;k++)
				b[i][j]=b[i][j]+u[i][k]*l[k][j];
		}
	}

    return 1;
}