/*  prog2c.m -- Cholesky Factorization
    Written Oct 10, 2007 by Eric Olson for Math/CS 466/666

    For more information see Step 3 in First Steps in Numerical
    Analysis, 2nd Edition, by Hosking, Joe, Joyce and Turner */

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

static FILE *fip;

typedef struct {
    double *d;
    int n;
} matrix;

matrix getmatrix(){
    int i,j,r;
    static matrix A={0,0};
    r=fscanf(fip,"%d",&A.n);
    if(r!=1) A.n=0;
    if(A.n==0) return A;
    if(A.d) free(A.d);
    A.d=(double *)malloc(sizeof(double)*A.n*A.n);
    for(i=0;i<A.n;i++) for(j=0;j<A.n;j++) {
        r=fscanf(fip,"%lg",&A.d[i*A.n+j]);
        if(r!=1) {
            A.n=0;
            return A;
        }
    }
    return A;
}

void putmatrix(matrix A){
    int i,j;
    for(i=0;i<A.n;i++) {
        for(j=0;j<A.n;j++){
            printf(" %10g",A.d[i*A.n+j]);
        }
        printf("\n");
    }
}

int main(){
    matrix A;
    printf("prog2c.m -- Cholesky Factorization\n\n");
    fip=fopen("prog2c.dat","r");
    if(!fip){
        perror("main");
        exit(-1);
    }
    for(;;)
    {
        A=getmatrix();
        if(A.n==0) break;
        putmatrix(A);
        printf("\n");

        /* add code here */

    }
    fclose(fip);
    exit(0);
}