/**
   * Lényege, hogy 1 mátrixnak előállítjuk az inkomplett cholesky felbontását,
   * amelyet a prekondícionált konjugált mátrixnál használunk fel, mint prekondícionálási
   * mátrixot.
   * A mátrix szimmetrikus és pozitív definit kell legyen.
   * felteszem, hogy az, mert csak ilyenekkel dolgozunk.
   */
   /* Megoldás a könyvbeni pszeudo kóddal. */

void Inkomplett_Cholesky();

#include <stdio.h>
#include <math.h>
#include <stdlib.h>
/*#include <sys\timeb.h>*/

int main (void){
    /*struct timeb t_start, t_current;
    int t_diff_in_sec, t_diff_in_ms;
    ftime(&t_start);*/
    
    Inkomplett_Cholesky();
    
    /*ftime(&t_current);
      t_diff_in_sec = (int)  (t_current.time - t_start.time);
      t_diff_in_ms = (int ) (t_current.millitm - t_start.millitm);

    printf("\nEltelt idő: %d ms",t_diff_in_sec*1000+t_diff_in_ms);*/
    return 0;
}
/*____________________________________________________________________________*/
void Inkomplett_Cholesky(){
    int n, i, j, k;
    n = i = j = k = 0;
    
    scanf("%d\n",&n);
    printf("n: %d\n", n);
    double (*A)[n] = NULL;
    A = calloc (n, sizeof (*A));

    /* A matrix beolvasasa */
    for ( i=0; i<n; i++ ){
        for ( j=0; j<n; j++ ){
            scanf("%lf", &A[i][j]);
        }
    }

    
    /* Az algoritmus */
    for ( k=0; k<n; k++){
      A[k][k] = sqrt( A[k][k] );
      for( i=k+1; i<n; i++){
        if ( fabs( A[i][k] ) > pow ( 10,-15 ) ){
          /*printf( "%lf / %lf ", A[i][k], A[k][k]);*/
          A[i][k] = A[i][k] / A[k][k];
          /*printf("A[%d][%d] = %lf\n", i, k, A[i][k]);*/
        }
      }
      for( j=k+1; j<n; j++ ){
        for( i=j; i<n; i++ ){
          if ( fabs( A[i][j] ) > pow ( 10,-15 ) ){
            /*printf( "%lf - (%lf * %lf) ", A[i][j], A[i][k], A[j][k] );
            printf("A[%d][%d] = %lf\n", i, j, A[i][j]);*/
            A[i][j] = A[i][j] - ( A[i][k] * A[j][k] );
            /*printf("A[%d][%d] = %lf\n", i, j, A[i][j]);*/
          }
        }
      }
    }
    
    for ( i=0; i<n; i++){
      for( j=i+1; j<n; j++ ){
        A[i][j] = 0.0;
      }
    }
    
    puts("");
    for ( i=0; i<n; i++){
      for( j=0; j<n; j++  ){
        printf( "%lf ", A[i][j] );
      }
      puts("");
    }
    
return;
}
