1 /* cminv.c CCMATH mathematics library source code.
2 *
3 * Copyright (C) 2000 Daniel A. Atkinson All rights reserved.
4 * This code may be redistributed under the terms of the GNU library
5 * public license (LGPL). ( See the lgpl.license file for details.)
6 * ------------------------------------------------------------------------
7 */
8 #include <stdlib.h>
9 #include "complex.h"
cminv(Cpx * a,int n)10 int cminv(Cpx *a,int n)
11 { int i,j,k,m,lc,*le; Cpx *ps,*p,*q,*pa,*pd;
12 Cpx z,h,*q0; double s,t,tq=0.,zr=1.e-15;
13 le=(int *)calloc(n,sizeof(int));
14 q0=(Cpx *)calloc(n,sizeof(Cpx));
15 pa=pd=a;
16 for(j=0; j<n ;++j,++pa,pd+=n+1){
17 if(j>0){
18 for(i=0,p=pa,q=q0; i<n ;++i,p+=n) *q++ = *p;
19 for(i=1; i<n ;++i){ lc=i<j?i:j;
20 z.re=z.im=0.;
21 for(k=0,p=pa+i*n-j,q=q0; k<lc ;++k,++q,++p){
22 z.re+=p->re*q->re-p->im*q->im;
23 z.im+=p->im*q->re+p->re*q->im;
24 }
25 q0[i].re-=z.re; q0[i].im-=z.im;
26 }
27 for(i=0,p=pa,q=q0; i<n ;++i,p+=n) *p= *q++;
28 }
29 s=fabs(pd->re)+fabs(pd->im); lc=j;
30 for(k=j+1,ps=pd; k<n ;++k){ ps+=n;
31 if((t=fabs(ps->re)+fabs(ps->im))>s){ s=t; lc=k;}
32 }
33 tq=tq>s?tq:s; if(s<zr*tq){ free(le-j); free(q0); return -1;}
34 *le++ =lc;
35 if(lc!=j){ p=a+n*j; q=a+n*lc;
36 for(k=0; k<n ;++k,++p,++q){ h= *p; *p= *q; *q=h;}
37 }
38 t=pd->re*pd->re+pd->im*pd->im;
39 h.re=pd->re/t; h.im= -(pd->im)/t;
40 for(k=j+1,ps=pd+n; k<n ;++k,ps+=n){
41 z.re=ps->re*h.re-ps->im*h.im;
42 z.im=ps->im*h.re+ps->re*h.im; *ps=z;
43 }
44 *pd=h;
45 }
46 for(j=1,pd=ps=a; j<n ;++j){
47 for(k=0,pd+=n+1,q= ++ps; k<j ;++k,q+=n){
48 z.re=q->re*pd->re-q->im*pd->im;
49 z.im=q->im*pd->re+q->re*pd->im; *q=z;
50 }
51 }
52 for(j=1,pa=a; j<n ;++j){ ++pa;
53 for(i=0,q=q0,p=pa; i<j ;++i,p+=n) *q++ = *p;
54 for(k=0; k<j ;++k){ h.re=h.im=0.;
55 for(i=k,p=pa+k*n+k-j,q=q0+k; i<j ;++i){
56 h.re-=p->re*q->re-p->im*q->im;
57 h.im-=p->im*q->re+p->re*q->im; ++p; ++q;
58 }
59 q0[k]=h;
60 }
61 for(i=0,q=q0,p=pa; i<j ;++i,p+=n) *p= *q++;
62 }
63 for(j=n-2,pd=pa=a+n*n-1; j>=0 ;--j){ --pa; pd-=n+1;
64 for(i=0,m=n-j-1,q=q0,p=pd+n; i<m ;++i,p+=n) *q++ = *p;
65 for(k=n-1,ps=pa; k>j ;--k,ps-=n){
66 z.re= -ps->re; z.im= -ps->im;
67 for(i=j+1,p=ps+1,q=q0; i<k ;++i,++p,++q){
68 z.re-=p->re*q->re-p->im*q->im;
69 z.im-=p->im*q->re+p->re*q->im;
70 }
71 q0[--m]=z;
72 }
73 for(i=0,m=n-j-1,q=q0,p=pd+n; i<m ;++i,p+=n) *p= *q++;
74 }
75 for(k=0,pa=a; k<n-1 ;++k,++pa){
76 for(i=0,q=q0,p=pa; i<n ;++i,p+=n) *q++ = *p;
77 for(j=0,ps=a; j<n ;++j,ps+=n){
78 if(j>k){ h.re=h.im=0.; p=ps+j; i=j;}
79 else{ h=q0[j]; p=ps+k+1; i=k+1;}
80 for(; i<n ;++i,++p){
81 h.re+=p->re*q0[i].re-p->im*q0[i].im;
82 h.im+=p->im*q0[i].re+p->re*q0[i].im;
83 }
84 q0[j]=h;
85 }
86 for(i=0,q=q0,p=pa; i<n ;++i,p+=n) *p= *q++;
87 }
88 for(j=n-2,le--; j>=0 ;--j){
89 for(k=0,p=a+j,q=a+ *(--le); k<n ;++k,p+=n,q+=n){
90 h= *p; *p= *q; *q=h;
91 }
92 }
93 free(le); free(q0);
94 return 0;
95 }
96