text: #include "trace.h"
  #include "common.h"
  #include <math.h>
  #include <omp.h>
  #include <stdlib.h>
  
  void panel(double *P, int k, info_type info){
  
    int m, ld, err, np, i;
  
    i   = k*info.B;
    np  = info.B;
    m   = info.N-i;
    ld  = info.N;
  #if defined(DBG)
    printf("%2d -- panel : %d\n",omp_get_thread_num(),k);
  #endif
    trace_event_start(PNL);
    dgetrf_(&m, &np, P+i, &ld, info.ipiv+i, &err);
    trace_event_stop(PNL);
  }
  
  
  void update(double *P, double *U, int k, int j, info_type info){
  
    int m, mu, nu, np, ld, err, i;
    char NoTran = 'N', Lower='L', Unit='U', Left='L';
    int IONE=1, IMONE=-1, IZERO=0;
    double DONE=(double)1.0, DMONE=-1.0, DZERO=0.0;
  
    i  = k*info.B;
    np = info.B;
    nu = info.B;
  
    m   = info.N-i;
    mu  = m-np;
    ld  = info.N;
    
  #if defined(DBG)
    printf("%2d -- update: %d %d\n",omp_get_thread_num(),k,j);
  #endif
    trace_event_start(UPD);
    dlaswp_(&nu, U+i, &ld, &IONE, &np, info.ipiv+i, &IONE);
    dtrsm_ (&Left, &Lower, &NoTran, &Unit, 
            &np, &nu, 
            &DONE, 
            P+i, &ld,
            U+i, &ld);
    dgemm_ (&NoTran, &NoTran,
            &mu, &nu, &np,
            &DMONE, 
            P+i+np, &ld,
            U+i, &ld,
            &DONE,
            U+i+np, &ld);
    trace_event_stop(UPD);
  
    return;
  }
  
  
  void backperm(Matrix A, info_type info){
  
    int i, j, ld, ipb, ipo;
    int IONE=1;
  
    ld  = info.N;
  
    trace_event_start(END);
    for(i=info.B; i<info.N; i+=info.B){
      for(j=i; j<info.N; j++)
        info.ipiv[j]+=info.B;
      ipo = i+1;
      ipb = i+info.B;
  #pragma omp parallel for
        for(j=0; j<i/info.B; j+=1)
          dlaswp_(&info.B, A[j], &ld, &ipo, &ipb, info.ipiv, &IONE);
    }
    trace_event_stop(END);
  
    return;
  }
  
  
  void checkres(double *A, double *Acpy, info_type info){
  
    int     ISEED[4] = {0,0,0,3};
    int     IONE=1;
    char    NoTran = 'N';
    double  DONE=1.0, DMONE=-1.0;
    double  *x, *b;
    int err;
    double  flops, nrm2, nrmf;
  
    x      = (double *)malloc(info.N*sizeof(double));
    b      = (double *)malloc(info.N*sizeof(double));
  
  
    dlarnv_(&IONE, ISEED, &info.N,  b);
    memcpy( x, b, info.N*sizeof(double) );
  
    dgetrs_(&NoTran, &(info.N), &IONE, A, &(info.N), info.ipiv, x, &(info.N), &err);
  
    dgemv_(&NoTran, &info.N, &info.N, &DONE, Acpy, &info.N, x, &IONE, &DMONE, b, &IONE);
    nrm2 = dnrm2_c(info.N, b, IONE);
    nrmf = dnrmf_c(info.N, info.N, Acpy, info.N);
    printf("Residual norm: %e\n",sqrt(info.N)*nrm2/nrmf);
  
    free(x);
    free(b);
  
    return;
  
  }
  