/* cmatinv.cpp -- inverse of complex matrix using Gauss-Jordan elimination.
 * (C) 2001, C. Bond. All rights reserved.
 */
#include <complex.h>
/* Bare-bones Gauss-Jordan algorithm. Returns '0' on success, '1' if
 * matrix is singular. Only simple pivoting used.
 *
 * Replaces input matrix with its inverse.
 */
int cmatinv(complex<double> **a,int n)
{
    complex<double> **inv,tmp,*t;
    int i,j,k,retval;

    retval = 0;
    inv = new complex<double> *[n];
    for (i=0;i<n;i++) {
        inv[i] = new complex<double> [n];
    }
// Initialize identity matrix
    for (i=0;i<n;i++) {
        for (j=0;j<n;j++) {
            inv[i][j] = 0.0;
        }
        inv[i][i] = 1.0;
    }

    for (k=0;k<n;k++) {
        tmp = a[k][k];
        if (abs(tmp)==0.0) {
// pivoting
            for (j=k+1;j<n;j++) {
                if((tmp = abs(a[j][k])) != 0.0) break;
            }
            if (tmp == 0.0) {
                retval = 1;
                goto _100;
            }
            t=a[j];                 // swap matrix rows...
            a[j]=a[k];
            a[k]=t;
            t=inv[j];               // ...and result matrix
            inv[j]=inv[k];
            inv[k]=t;
        }
        for (j=0;j<n;j++) {
            if (j>k) a[k][j] /= tmp;    // Don't bother with previous entries
            inv[k][j] /= tmp;
        }
        for (i=0;i<n;i++) {             // Loop over rows
            if (i == k) continue;
            tmp = a[i][k];
            for (j=0;j<n;j++) {
                if (j>k) a[i][j] -= a[k][j]*tmp;
                inv[i][j] -= inv[k][j]*tmp;
            }
        }
    }

// Copy inverse to source matrix
    for (i=0;i<n;i++) {
        for (j=0;j<n;j++) {
            a[i][j] = inv[i][j];
        }
    }
_100:
    for (i=0;i<n;i++) {
        delete [] inv[i];
    }
    delete [] inv;
    return retval;
}