C 求解反矩陣的程式

Question

最近做論文遇到要求解反矩陣
因為我論文的程式是用C 撰寫，不像 Matlab可以呼叫function
請問誰也C code??

Answer 1

這些是我兩年前寫的程式，沒有最佳化，沒有寫 permutation 的那部份。
用法請參考 main( )。
除了 main 裡要的以外，其它需要的副程式應該都有了。
要是真的欠了什麼，再 email 告訴我。
/*Inverse of a matrix
byJacob Lee7/30/2005for Master thesis
ma_inv(**,) returns a pointer to the calculated matrix or
NULL if matrix is singular or
memory is not enough
This version does not do permutation!!
It may misjudge a non-singular matrix to be singularand
the value may less accurate!!
The caller must free both (*t) and (t)
include dim.cor dim.oand
matcopy.cor matcopy.oand
mateye.cor mateye.o and
matpermu.cor matpermu.oin its program or project.
See main() for the usage, especially, the order of free()!*/

#include
#include
#include //printf();
#include
#include"matrix.h"

MA_TYP**ma_inv(MA_TYP **A, int s)
{MA_TYPc, **I, **L, **U, **T, max, *tmp;//Cache
inti, j, p, t, *pp;
if (U= ma_copy(A, s, s))//if (I != NULL)
{if (L= ma_eye(s))
{//if (pp= ma_permu(s))
{for (i=0; i {/*max= fabs(U[i][i]);//Choice max
p= i;
for (j=i 1; j if (max < fabs(U[j][i]))
max= fabs(U[j][i]),p= j;
if (p > i)// need swap?
{tmp= (MA_TYP *) malloc(t= sizeof(*tmp) * (s - i));
memcpy(tmp, &U[i][i],t);//Swap U
memcpy(&U[i][i], &U[p][i],t);
memcpy(&U[p][i], tmp,t);
free(tmp);
tmp= (MA_TYP *) malloc(t= sizeof(*tmp) * i);
memcpy(tmp,&L[i][0],t);//Swap L
memcpy(&L[i][0],&L[p][0],t);
memcpy(&L[p][0],tmp,t);
free(tmp);
t = pp[i],pp[i] = pp[p],pp[p] = t;
}*/

2007-03-08 09:39:41 補充：
// make the lower part 0.
// max is now useless, use it as tmp
if (c = U[i][i])//non-singular?
{ for (j=i+1; j<s; j++)
{ max= U[j][i] / c;
for (t=i; t<s; t++)//U[i][0] ~ U[0][i-1] is 0.
U[j][t]-= U[i][t] * max;
for (t=0; t<s; t++)
L[j][t]-= L[i][t] * max;
}

2007-03-08 09:40:09 補充：
//make the upperer part 0?
for (j=0; j<i; j++)
{max= U[j][i] / c;
//U[j][i]= 0.;//useless!!!
for (t=i+1; t<s; t++)//U[i][0] ~ U[0][i-1] is 0.
U[j][t]-= U[i][t] * max;
for (t=0;t<s; t++)
L[j][t]-= L[i][t] * max;
}

2007-03-08 09:45:43 補充：
// normalize ith row. Why I can not do these later?
// U[i][i]= 1.; // useless!!!
for (t=i+1; t<s; t++) // U[i][0] ~ U[0][i-1] is 0.
U[i][t]/= c;
for (t=0;t<s; t++)
L[i][t]/= c;
}
else // singular!
{ free(*L); free(L); L = 0; break; }
} }

2007-03-08 09:49:02 補充：
// else // Because no permu this version
// { free(*L); free(L); L = 0; } //Because no Permu this version
}
// else // I think this is a bug!
// { free(*U); free(U);} // I think this is a bug!
free(*U);free(U);
}
returnL;
}

2007-03-08 09:54:34 補充：
int size_of[] = { sizeof(char), sizeof(unsigned char), sizeof(int), sizeof(unsigned),
sizeof(long), sizeof(unsigned long), sizeof(float), sizeof(double), sizeof(long double) };

enum {D_CHAR, D_U_CHAR, D_INT, D_UNSIGNED, D_LONG, D_U_LONG, D_FLOAT, D_DOUBLE, D_L_DOUBLE}

2007-03-08 09:56:20 補充：
void **dim2(void ***A, int h, int w, int type)
{ int i, n;
char *buf, **idx;

*A = 0; // Assume allocation fail
if ( idx = (char **) malloc(sizeof(char*) * h) )//if (!NULL)
if ( buf = (char *) malloc(size_of[type] * h * w) )
{ w *= size_of[type];

2007-03-08 09:56:51 補充：
for (n=i=0; i<h; i++, n+=w)
idx[i] = &buf[n];

*A = (void**) idx;
}
else free(idx); // Can not alloc the 2nd, free the 1st

return *A;
}

2007-03-08 10:01:36 補充：
MA_TYP**ma_copy(MA_TYP **A, int h, int w)
{ MA_TYP **C;

if (dim2((void ***)&C, h, w, DI_TYP)) // if (!NULL)
memcpy(*C, *A, h * w * sizeof(**A));

returnC;
}

你要另外寫個
typedef double MA_TYP; // Change to float if not enough memory

除了 // 外，我用的是 C，所以沒有 function overload

2007-03-08 10:04:27 補充：
double **ma_eye(int s) // size
{ MA_TYP **I;
int i;

if (dim2((void ***)&I, s, s, DI_TYP)) // if (!NULL)
{ bzero(*I, sizeof(**I) * s * s); // VC 沒有 bzero，你就用 for loop 把它清為 0
for (i=0; i<s; i++)
I[i][i] = 1.;
}

return I;
}

2007-03-08 10:08:01 補充：
int main(int argc, char **argv)
{MA_TYP **A, **B, **C, **I, **in, norm;
double CLK, d; // #(Double) to count
clock_t c, c0, c1, t, t0, t1; // clock, time
int i, iter, s, n; // Number of singular
char nm[20];

iter=argc==2?atoi(argv[1]):9;

srand48(time(0));
CLK=sysconf(_SC_CLK_TCK);

2007-03-08 10:10:54 補充：
t0=t1=clock();
c0=c1=times(0);

printf("Begin %d iterations\n", iter);
for (d = n = i = 0; i<iter; i++)
{s= lrand48() % 600 + 2;
//s= 1400+i; // used for basic testing
d += s * s;

//printf("Iter %4d / %d Size: [%4d]^2 \n", i, iter, s);//for basic test

2007-03-08 10:15:39 補充：
if (A = ma_rand(s, s, 20., 10.)) // if (A != NULL)
{//ma_trnc(&A, s, s); // Truncate to int for human easy to check
// free(*A); free(A); // for gaussian test
// A = gaussian(s, 1., 0); // for gaussian test

// Insert specified matrix, sampled at the bottom of this file, here if need

2007-03-08 10:18:27 補充：
//Be sure to change the value of s 10 lines above according to the matrix

// if (s<32) ma_show(A, s, s);// for showing if you want

if (in= ma_inv(A, s))
{if (I= ma_eye(s))
{if (B= ma_iner(A, in, s, s, s))
{if (C= ma_minu(B, I, s, s))
{norm= ma_norm(C, s, s);
if (norm > 1+0*EPSILON*2048)

2007-03-08 10:22:58 補充：
{printf("\nThe original %d * %d matrix, A, is\n", s, s);

ma_show(A, s, s);
printf("\nThe inverse matrix of A is\n");
ma_show(in, s, s);
printf("The 2 norm of A * inv(A) - I(A) is %lg%c\n", norm, 7);
}
free(*C); free(C);
}
free(*B); free(B);
}
if (B = ma_iner(in, A, s, s, s))

2007-03-08 10:28:50 補充：
if(C=ma_minu(B,I,s,s))
{norm=ma_norm(C,s,s);
if(norm>1+0*EPSILON*2048)
{printf("\nThe original %d * %d matrix, A, is\n",s,s);
ma_show(A,s,s);
printf("\nThe inverse matrix of A is\n");
ma_show(in,s,s);
printf("The 2 norm of I(A) - A * inv(A) is %lg%c\n",norm,7);
}

2007-03-08 10:31:48 補充：
free(*C); free(C);
}
free(*B); free(B);
}
free(*I); free(I);
}
free(*in); free(in);
}
else
{//printf("Matrix A may be singular\n");
//ma_show(A, s, s);
n++;
} }
t= clock();
c= times(0);
if (in) sprintf(nm, "%lf", norm);

2007-03-08 10:33:33 補充：
printf("Iter %4d / %d Size: [%4d]^2 time used: %6.2lf / %7.2lf seconds %s is %s.\n",
i, iter, s, (t-t1) / (double)CLOCKS_PER_SEC, (c-c1)/CLK, in ? "norm" : "A",in ? nm : "singular");
t1 = t, c1 = c;
}

2007-03-08 10:34:14 補充：
printf("Complete %d iterations in %.2lf / %.2lf seconds, sum(matrix size) = %.0lf. %d of them %s singular%s\n",
iter, (t-t0) / (double) CLOCKS_PER_SEC, (c-c0)/CLK,d, n, n<2 ? "is" : "are", n<2 ? "." : "s.");

return 0;
}