#include "PLA.h"

int FLA_Chol_unb_var3( PLA_Obj A )
{
  PLA_Obj ATL=NULL,   ATR=NULL,      A00=NULL,  a01=NULL,     A02=NULL, 
          ABL=NULL,   ABR=NULL,      a10t=NULL, alpha11=NULL, a12t=NULL,
                                     A20=NULL,  a21=NULL,     A22=NULL;

  PLA_Obj MINUS_ONE=NULL, ZERO=NULL, ONE=NULL;

  PLA_Create_constants_conf_to( A, &MINUS_ONE, &ZERO, &ONE );

  PLA_Part_2x2( A,    &ATL, &ATR,
                      &ABL, &ABR,     0, 0, PLA_TL );

  while ( PLA_Obj_length( ATL ) < PLA_Obj_length( A ) ){

    PLA_Repart_2x2_to_3x3( ATL, /**/ ATR,       &A00,  /**/ &a01,     &A02,
                        /* ************* */   /* ************************** */
                                                &a10t, /**/ &alpha11, &a12t,
                           ABL, /**/ ABR,       &A20,  /**/ &a21,     &A22,
                           1, 1, PLA_BR );

    /*------------------------------------------------------------*/

    PLA_Sqrt( alpha11 );

    PLA_Inv_scal( alpha11, a21 );

    PLA_Syr( PLA_LOWER_TRIANGULAR, MINUS_ONE, a21, A22 );

    /*------------------------------------------------------------*/

    PLA_Cont_with_3x3_to_2x2( &ATL, /**/ &ATR,       A00,  a01,     /**/ A02,
                                                     a10t, alpha11, /**/ a12t,
                            /* ************** */  /* ************************ */
                              &ABL, /**/ &ABR,       A20,  a21,     /**/ A22,
                              PLA_TL );

  }

  PLA_Obj_free( &ATL ); PLA_Obj_free( &ATR );
  PLA_Obj_free( &ABL ); PLA_Obj_free( &ABR );
  PLA_Obj_free( &A00 );  PLA_Obj_free( &a01 );     PLA_Obj_free( &A02 );
  PLA_Obj_free( &a10t ); PLA_Obj_free( &alpha11 ); PLA_Obj_free( &a12t );
  PLA_Obj_free( &A20 );  PLA_Obj_free( &a21 );     PLA_Obj_free( &A22 );

  PLA_Obj_free( &MINUS_ONE ); PLA_Obj_free( &ZERO ); PLA_Obj_free( &ONE );

  return PLA_SUCCESS;
}