text:
/*
 * Copyright (c) 2011-2018 The University of Tennessee and The University
 *                         of Tennessee Research Foundation.  All rights
 *                         reserved.
 *
 * @precisions normal z -> s d c
 *
 */
#include "dplasma_cores.h"
#include "dplasma_zcores.h"

#if defined(PARSEC_HAVE_STRING_H)
#include <string.h>
#endif  /* defined(PARSEC_HAVE_STRING_H) */
#if defined(PARSEC_HAVE_STDARG_H)
#include <stdarg.h>
#endif  /* defined(PARSEC_HAVE_STDARG_H) */
#include <stdio.h>
#ifdef PARSEC_HAVE_LIMITS_H
#include <limits.h>
#endif
#include <stdlib.h>

#include <cblas.h>
#include <core_blas.h>

#define max(a, b) ((a) > (b) ? (a) : (b))
#define min(a, b) ((a) < (b) ? (a) : (b))
#if defined(ADD_)
#define  DLARFG      zlarfg_
#else
#define  DLARFG      zlarfg
#endif
extern void DLARFG(int *N, parsec_complex64_t *ALPHA, parsec_complex64_t *X, int *INCX, parsec_complex64_t *TAU);

//static void band_to_trd_vmpi1(int N, int NB, parsec_complex64_t *A, int LDA);
//static void band_to_trd_vmpi2(int N, int NB, parsec_complex64_t *A, int LDA);
//static void band_to_trd_v8seq(int N, int NB, parsec_complex64_t *A, int LDA, int INgrsiz, int INthgrsiz);
//static int  TRD_seqgralgtype(int N, int NB, parsec_complex64_t *A, int LDA, parsec_complex64_t *C, parsec_complex64_t *S, int i, int j, int m, int grsiz, int BAND);
int  blgchase_ztrdv1(int NT, int N, int NB, parsec_complex64_t *A, parsec_complex64_t *V, parsec_complex64_t *TAU, int sweep, int id, int blktile);
int  blgchase_ztrdv2(int NT, int N, int NB, parsec_complex64_t *A1, parsec_complex64_t *A2, parsec_complex64_t *V1, parsec_complex64_t *TAU1, parsec_complex64_t *V2, parsec_complex64_t *TAU2, int sweep, int id, int blktile);

int CORE_zlarfx2(int side, int N,
                parsec_complex64_t V,
                parsec_complex64_t TAU,
                parsec_complex64_t *C1, int LDC1,
                parsec_complex64_t *C2, int LDC2);

int CORE_zlarfx2c(int uplo,
                parsec_complex64_t V,
                parsec_complex64_t TAU,
                parsec_complex64_t *C1,
                parsec_complex64_t *C2,
                parsec_complex64_t *C3);

static void CORE_zhbtelr(int N, int NB, parsec_complex64_t *A1, int LDA1, parsec_complex64_t *A2, int LDA2, parsec_complex64_t *V1, parsec_complex64_t *TAU1, int st, int ed);
static void CORE_zhbtrce(int N, int NB, parsec_complex64_t *A1, int LDA1, parsec_complex64_t *A2, int LDA2, parsec_complex64_t *V1, parsec_complex64_t *TAU1, parsec_complex64_t *V2, parsec_complex64_t *TAU2, int st, int ed, int edglob);
static void CORE_zhbtlrx(int N, int NB, parsec_complex64_t *A1, int LDA1, parsec_complex64_t *A2, int LDA2, parsec_complex64_t *V1, parsec_complex64_t *TAU1, int st, int ed);

static void DLARFX_C(char side, int N, parsec_complex64_t V, parsec_complex64_t TAU, parsec_complex64_t *C, int LDC);
static void TRD_type1bHL(int N, int NB, parsec_complex64_t *A, int LDA, parsec_complex64_t *V, parsec_complex64_t *TAU, int st, int ed);
static void TRD_type2bHL(int N, int NB, parsec_complex64_t *A, int LDA, parsec_complex64_t *V, parsec_complex64_t *TAU, int st, int ed);
static void TRD_type3bHL(int N, int NB, parsec_complex64_t *A, int LDA, parsec_complex64_t *V, parsec_complex64_t *TAU, int st, int ed);

int CORE_zlarfx2(int side, int N,
                parsec_complex64_t V,
                parsec_complex64_t TAU,
                parsec_complex64_t *C1, int LDC1,
                parsec_complex64_t *C2, int LDC2)
{
    static parsec_complex64_t zzero = 0.0;
    int    J;
    parsec_complex64_t V2, T2, SUM;

    /* Quick return */
    /*
    if (N == 0)
        return PLASMA_SUCCESS;
    */
    if (TAU == zzero)
        return 0;

    /*
     * Special code for 2 x 2 Householder where V1 = I
    */
    if(side==PlasmaLeft){
        V2 = conj(V);
        T2 = TAU*conj(V2);
        for (J = 0; J < N ; J++){
           SUM         = C1[J*LDC1]   + V2*C2[J*LDC2];
           C1[J*LDC1]  = C1[J*LDC1]   - SUM*TAU;
           C2[J*LDC2]  = C2[J*LDC2]   - SUM*T2;
        }
    }else if(side==PlasmaRight){
        V2 = V;
        T2 = TAU*conj(V2);
        for (J = 0; J < N ; J++){
           SUM     = C1[J]   + V2*C2[J];
           C1[J]   = C1[J]   - SUM*TAU;
           C2[J]   = C2[J]   - SUM*T2;
        }
    }

    return 0;
}

/***************************************************************************//**
 *
 **/
int CORE_zlarfx2c(int uplo,
                parsec_complex64_t V,
                parsec_complex64_t TAU,
                parsec_complex64_t *C1,
                parsec_complex64_t *C2,
                parsec_complex64_t *C3)
{
    static parsec_complex64_t zzero = 0.0;
    parsec_complex64_t T2, SUM, TEMP;

    /* Quick return */
    if (TAU == zzero)
        return 0;

   /*
    *        Special code for a diagonal block  C1
    *                                           C2  C3
    */
 if(uplo==PlasmaLower){   //do the corner Left then Right (used for the lower case tridiag)
      // L and R for the 2x2 corner
      //            C(N-1, N-1)  C(N-1,N)        C1  TEMP
      //            C(N  , N-1)  C(N  ,N)        C2  C3
      // For Left : use conj(TAU) and V.
      // For Right: nothing, keep TAU and V.
      // Left 1 ==> C1
      //            C2
      TEMP    = conj(C2[0]); // copy C2 here before modifying it.
      T2      = conj(TAU)*V;
      SUM     = C1[0]   + conj(V)*C2[0];
      C1[0]   = C1[0]   - SUM*conj(TAU);
      C2[0]   = C2[0]   - SUM*T2;
      // Left 2 ==> TEMP
      //            C3
      SUM     = TEMP    + conj(V)*C3[0];
      TEMP    = TEMP    - SUM*conj(TAU);
      C3[0]   = C3[0]   - SUM*T2;
      // Right 1 ==>  C1 TEMP.  NB: no need to compute corner (2,2)=TEMP
      T2      = TAU*conj(V);
      SUM     = C1[0]   + V*TEMP;
      C1[0]   = C1[0]   - SUM*TAU;
      // Right 2 ==> C2 C3
      SUM     = C2[0]   + V*C3[0];
      C2[0]   = C2[0]   - SUM*TAU;
      C3[0]   = C3[0]   - SUM*T2;
 }else if(uplo==PlasmaUpper){ // do the corner Right then Left (used for the upper case tridiag)
      //            C(N-1, N-1)  C(N-1,N)        C1    C2
      //            C(N  , N-1)  C(N  ,N)        TEMP  C3
      // For Left : use TAU       and conj(V).
      // For Right: use conj(TAU) and conj(V).
      // Right 1 ==> C1 C2
      V       = conj(V);
      TEMP    = conj(C2[0]); // copy C2 here before modifying it.
      T2      = conj(TAU)*conj(V);
      SUM     = C1[0]   + V*C2[0];
      C1[0]   = C1[0]   - SUM*conj(TAU);
      C2[0]   = C2[0]   - SUM*T2;
      // Right 2 ==> TEMP C3
      SUM     = TEMP    + V*C3[0];
      TEMP    = TEMP    - SUM*conj(TAU);
      C3[0]   = C3[0]   - SUM*T2;
      // Left 1 ==> C1
      //            TEMP. NB: no need to compute corner (2,1)=TEMP
      T2      = TAU*V;
      SUM     = C1[0]   + conj(V)*TEMP;
      C1[0]   = C1[0]   - SUM*TAU;
      // Left 2 ==> C2
      //            C3
      SUM     = C2[0]   + conj(V)*C3[0];
      C2[0]   = C2[0]   - SUM*TAU;
      C3[0]   = C3[0]   - SUM*T2;
 }

 return 0;
}

///////////////////////////////////////////////////////////
//                  DLARFX en C
///////////////////////////////////////////////////////////
static void DLARFX_C(char side, int N, parsec_complex64_t V, parsec_complex64_t TAU, parsec_complex64_t *C, int LDC)
{
 parsec_complex64_t T2, SUM, TEMP;
 int    J, pt;

/*
 *        Special code for 2 x 2 Householder
 */

 T2 = TAU*V;
 if(side=='L'){
     for (J = 0; J < N ; J++){
        pt      = LDC*J;
        SUM     = C[pt]   + V*C[pt+1];
        C[pt]   = C[pt]   - SUM*TAU;
        C[pt+1] = C[pt+1] - SUM*T2;
     }
 }else if(side=='R'){
     for (J = 0; J < N ; J++){
        pt      = J+LDC;
        SUM     = C[J]   + V*C[pt];
        C[J]    = C[J]   - SUM*TAU;
        C[pt]   = C[pt]  - SUM*T2;
     }
 }else if(side=='B'){
     TEMP    = C[ LDC * (N-2) +1];
     for (J = 0; J < N-1 ; J++){
        pt      = LDC*J;
        SUM     = C[pt]   + V*C[pt+1];
        C[pt]   = C[pt]   - SUM*TAU;
        C[pt+1] = C[pt+1] - SUM*T2;
     }
     // L and R for the 2x2 corner
     //            C(1, N-1)  C(1,N)        C(1, N-1)  TEMP
     //            C(2, N-1)  C(2,N)        C(2, N-1)  C(2,N)
     // Left 1 ==> C(1,N-1) C(2,N-1) deja fait dans la boucle
     //pt      = LDC * (N-2);
     //TEMP    = C[pt+1];
     //SUM     = C[pt]   + V*C[pt+1];
     //C[pt]   = C[pt]   - SUM*TAU;
     //C[pt+1] = C[pt+1] - SUM*T2;
     // Left 2 ==> TEMP C(2,N)
     pt      = LDC * (N-1) +1;
     SUM     = TEMP    + V*C[pt];
     TEMP    = TEMP    - SUM*TAU;
     C[pt]   = C[pt]   - SUM*T2;
     // Right 1 ==>  C(1,N-1) TEMP.  NB: no need to compute corner (2,2)
     J       = LDC * (N-2);
     SUM     = C[J]   + V*TEMP;
     C[J]    = C[J]   - SUM*TAU;
     // Right 2 ==> C(2,N-1) C(2,N)
     J      = LDC * (N-2) + 1;
     pt     = LDC * (N-1) + 1;
     SUM    = C[J]   + V*C[pt];
     C[J]   = C[J]   - SUM*TAU;
     C[pt]  = C[pt]  - SUM*T2;
 }
}
///////////////////////////////////////////////////////////

///////////////////////////////////////////////////////////
#define A1(m,n)   &(A1[((m)-(n)) + LDA1*(n)])
#define A2(m,n)   &(A2[((m)-(n)) + LDA2*((n)-NB)])
#define V1(m)     &(V1[m-st])
#define TAU1(m)   &(TAU1[m-st])
#define V2(m)     &(V2[m-st])
#define TAU2(m)   &(TAU2[m-st])
///////////////////////////////////////////////////////////
//                  TYPE 1-BAND Householder
///////////////////////////////////////////////////////////
static void CORE_zhbtelr(int N, int NB, parsec_complex64_t *A1, int LDA1, parsec_complex64_t *A2, int LDA2, parsec_complex64_t *V1, parsec_complex64_t *TAU1, int st, int ed) {
  int    J1, J2, KDM1, LDX;
  int    len, len1, len2, t1ed, t2st;
  int    i, IONE, ITWO;
  IONE=1;
  ITWO=2;
  (void)N;

  KDM1 = NB-1;
  LDX = LDA2-1;
  /* **********************************************************************************************
   *   Annihiliate, then LEFT:
   * ***********************************************************************************************/
  for (i = ed; i >= st+1 ; i--){
     /* generate Householder to annihilate a(i+k-1,i) within the band */
     *V1(i)          = *A1(i, (st-1));
     *A1(i, (st-1))  = 0.0;
     DLARFG( &ITWO, A1((i-1),(st-1)), V1(i), &IONE, TAU1(i) );

     J1   = st;
     J2   = i-2;
     t1ed = min(J2,KDM1);
     t2st = max(J1, NB);
     len1 = t1ed - J1 +1;
     len2 = J2 - t2st +1;
     // printf("Type 1L st %d   ed %d    i %d    J1 %d   J2 %d  t1ed %d  t2st %d len1 %d len2 %d \n", st, ed, i, J1,J2,t1ed,t2st,len1,len2);
     /* apply reflector from the left (horizontal row) and from the right for only the diagonal 2x2.*/
     if(len2>=0){
        /* part of the left(if len2>0) and the corner are on tile T2 */
        if(len2>0) CORE_zlarfx2(PlasmaLeft, len2 , *V1(i), conj(*TAU1(i)), A2((i-1), t2st), LDX, A2(i, t2st), LDX);
        CORE_zlarfx2c(PlasmaLower, *V1(i), *TAU1(i), A2(i-1,i-1), A2(i,i-1), A2(i,i));
        if(len1>0) CORE_zlarfx2(PlasmaLeft, len1 , *V1(i), conj(*TAU1(i)), A1(i-1, J1), LDX, A1(i, J1), LDX);
     }else if(len2==-1){
        /* the left is on tile T1, and only A(i,i) of the corner is on tile T2 */
        CORE_zlarfx2c(PlasmaLower, *V1(i), *TAU1(i), A1(i-1,i-1), A1(i,i-1), A2(i,i));
        if(len1>0) CORE_zlarfx2(PlasmaLeft, len1 , *V1(i), conj(*TAU1(i)), A1(i-1, J1), LDX, A1(i, J1), LDX);
     }else{
        /* the left and the corner are on tile T1, nothing on tile T2 */
        CORE_zlarfx2c(PlasmaLower, *V1(i), *TAU1(i), A1(i-1,i-1), A1(i,i-1), A1(i,i)) ;
        if(len1>0) CORE_zlarfx2(PlasmaLeft, len1 , *V1(i), conj(*TAU1(i)), A1((i-1), J1), LDX, A1(i, J1), LDX);
     }
  }
  /* **********************************************************************************************
   *   APPLY RIGHT ON THE REMAINING ELEMENT OF KERNEL 1
   * ***********************************************************************************************/
  for (i = ed; i >= st+1 ; i--){
     J1    = i+1;
     J2    = ed;
     len   = J2-J1+1;
     if(len>0){
        if(i>NB)
           /* both column (i-1) and i are on tile T2 */
           CORE_zlarfx2(PlasmaRight, len, *V1(i), *TAU1(i), A2(J1, i-1), LDX, A2(J1, i), LDX);
        else if(i==NB)
           /* column (i-1) is on tile T1 while column i is on tile T2 */
           CORE_zlarfx2(PlasmaRight, len, *V1(i), *TAU1(i), A1(J1, i-1), LDX, A2(J1, i), LDX);
        else
           /* both column (i-1) and i are on tile T1 */
           CORE_zlarfx2(PlasmaRight, len, *V1(i), *TAU1(i), A1(J1, i-1), LDX, A1(J1, i), LDX);
     }
  }
}
///////////////////////////////////////////////////////////

///////////////////////////////////////////////////////////
//                  TYPE 2-BAND Householder
///////////////////////////////////////////////////////////
static void CORE_zhbtrce(int N, int NB, parsec_complex64_t *A1, int LDA1, parsec_complex64_t *A2, int LDA2, parsec_complex64_t *V1, parsec_complex64_t *TAU1, parsec_complex64_t *V2, parsec_complex64_t *TAU2, int st, int ed, int edglob) {
  int    J1, J2, J3, KDM1, LDX, pt;
  int    len, len1, len2, t1ed, t2st, iglob;
  int    i, IONE, ITWO;
  parsec_complex64_t V,T,SUM;
  IONE=1;
  ITWO=2;

  iglob = edglob+1;
  LDX   = LDA1-1;
  KDM1  = NB-1;
  /* **********************************************************************************************
   *   Right:
   * ***********************************************************************************************/
  for (i = ed; i >= st+1 ; i--){
     /* apply Householder from the right. and create newnnz outside the band if J3 < N */
     iglob = iglob -1;
     J1  = ed+1;
     if((iglob+NB)>= N){
        J2 = i +(N-iglob-1);
        J3 = J2;
     }else{
        J2 = i + KDM1;
        J3 = J2+1;
     }
     len   = J2-J1+1;
     /* printf("Type 2R st %d   ed %d    i %d   J1 %d   J2 %d  len %d iglob %d  \n",st,ed,i,J1,J2,len,iglob);*/

     if(len>0){
        if(i>NB){
           /* both column (i-1) and i are on tile T2 */
           CORE_zlarfx2(PlasmaRight, len, *V1(i), *TAU1(i), A2(J1, i-1), LDX, A2(J1, i), LDX);
           /* if nonzero element need to be created outside the band (if index < N) then create and eliminate it. */
           if(J3>J2){
              /* new nnz at TEMP=V2(i)  */
              V          = *V1(i);
              T          = *TAU1(i) * conj(V);
              SUM        =  V * (*A2(J3, i));
              *V2(i)     = -SUM * (*TAU1(i));
              *A2(J3, i) = *A2(J3, i) - SUM * T;
              /* generate Householder to annihilate a(j+kd,j-1) within the band */
              DLARFG( &ITWO, A2(J2,i-1), V2(i), &IONE, TAU2(i) );
           }
        }else if(i==NB){
           /* column (i-1) is on tile T1 while column i is on tile T2 */
           CORE_zlarfx2(PlasmaRight, len, *V1(i), *TAU1(i), A1(J1, i-1), LDX, A2(J1, i), LDX);
           /* if nonzero element need to be created outside the band (if index < N) then create and eliminate it. */
           if(J3>J2){
              /* new nnz at TEMP=V2(i)  */
              V          = *V1(i);
              T          = *TAU1(i) * conj(V);
              SUM        =  V * (*A2(J3, i));
              *V2(i)     = -SUM * (*TAU1(i));
              *A2(J3, i) = *A2(J3, i) - SUM * T;
              /* generate Householder to annihilate a(j+kd,j-1) within the band */
              DLARFG( &ITWO, A1(J2, i-1), V2(i), &IONE, TAU2(i) );
           }
        }else{
           /* both column (i-1) and i are on tile T1 */
           CORE_zlarfx2(PlasmaRight, len, *V1(i), *TAU1(i), A1(J1, i-1), LDX, A1(J1, i), LDX);
           /* if nonzero element need to be created outside the band (if index < N) then create and eliminate it. */
           if(J3>J2){
              /* new nnz at TEMP=V2(i)  */
              V              = *V1(i);
              T              = *TAU1(i) * conj(V);
              SUM            =  V * (*A1(J3, i));
              *V2(i)         = -SUM * (*TAU1(i));
              *A1(J3, i)    = *A1(J3, i) - SUM * T;
              /* generate Householder to annihilate a(j+kd,j-1) within the band */
              DLARFG( &ITWO, A1(J2, i-1), V2(i), &IONE, TAU2(i) );
           }
        }
     }
  }
              // if(id==1) return;

  /* **********************************************************************************************
   *   APPLY LEFT ON THE REMAINING ELEMENT OF KERNEL 1
   * ***********************************************************************************************/
  iglob = edglob+1;
  for (i = ed; i >= st+1 ; i--){
     iglob = iglob -1;
     if((iglob+NB)< N){ /* mean that J3>J2 and so a nnz has been created and so a left is required. */
        J1   = i;
        J2   = ed;
        t1ed = min(J2,KDM1);
        t2st = max(J1, NB);
        len1 = t1ed - J1 +1;
        len2 = J2 - t2st +1;
        pt   = i + KDM1; /* pt correspond to the J2 position of the corresponding right done above */
        //printf("Type 2L st %d   ed %d    i %d    J1 %d   J2 %d  t1ed %d  t2st %d len1 %d len2 %d \n", st, ed, i, J1,J2,t1ed,t2st,len1,len2);

        /* apply reflector from the left (horizontal row) and from the right for only the diagonal 2x2.*/
        if(len2>0){
           /* part of the left(if len2>0) and the corner are on tile T2 */
           CORE_zlarfx2(PlasmaLeft, len2 , *V2(i), conj(*TAU2(i)), A2(pt, t2st), LDX, A2(pt+1, t2st), LDX);
           if(len1>0) CORE_zlarfx2(PlasmaLeft, len1 , *V2(i), conj(*TAU2(i)), A1(pt, J1), LDX, A1(pt+1, J1), LDX);
        }else if(len1>0){
           /* the left and the corner are on tile T1, nothing on tile T2 */
           CORE_zlarfx2(PlasmaLeft, len1 , *V2(i), conj(*TAU2(i)), A1(pt, J1), LDX, A1(pt+1, J1), LDX);
        }
     }
  }
}
///////////////////////////////////////////////////////////

///////////////////////////////////////////////////////////
//                  TYPE 1-BAND Householder
///////////////////////////////////////////////////////////
static void CORE_zhbtlrx(int N, int NB, parsec_complex64_t *A1, int LDA1, parsec_complex64_t *A2, int LDA2, parsec_complex64_t *V1, parsec_complex64_t *TAU1, int st, int ed) {
  int    J1, J2, KDM1, LDX;
  int    len, len1, len2, t1ed, t2st;
  int    i;

  (void)N;
  KDM1 = NB-1;
  LDX = LDA2-1;
  /* **********************************************************************************************
   *   Annihiliate, then LEFT:
   * ***********************************************************************************************/
  for (i = ed; i >= st+1 ; i--){
     J1   = st;
     J2   = i-2;
     t1ed = min(J2,KDM1);
     t2st = max(J1, NB);
     len1 = t1ed - J1 +1;
     len2 = J2 - t2st +1;
     //printf("Type 3L st %d   ed %d    i %d    J1 %d   J2 %d  t1ed %d  t2st %d len1 %d len2 %d \n", st, ed, i, J1,J2,t1ed,t2st,len1,len2);
     /* apply reflector from the left (horizontal row) and from the right for only the diagonal 2x2.*/
     if(len2>=0){
        /* part of the left(if len2>0) and the corner are on tile T2 */
        if(len2>0) CORE_zlarfx2(PlasmaLeft, len2 , *V1(i), conj(*TAU1(i)), A2((i-1), t2st), LDX, A2(i, t2st), LDX);
        CORE_zlarfx2c(PlasmaLower, *V1(i), *TAU1(i), A2(i-1,i-1), A2(i,i-1), A2(i,i));
        if(len1>0) CORE_zlarfx2(PlasmaLeft, len1 , *V1(i), conj(*TAU1(i)), A1(i-1, J1), LDX, A1(i, J1), LDX);
     }else if(len2==-1){
        /* the left is on tile T1, and only A(i,i) of the corner is on tile T2 */
        CORE_zlarfx2c(PlasmaLower, *V1(i), *TAU1(i), A1(i-1,i-1), A1(i,i-1), A2(i,i));
        if(len1>0) CORE_zlarfx2(PlasmaLeft, len1 , *V1(i), conj(*TAU1(i)), A1(i-1, J1), LDX, A1(i, J1), LDX);
     }else{
        /* the left and the corner are on tile T1, nothing on tile T2 */
        CORE_zlarfx2c(PlasmaLower, *V1(i), *TAU1(i), A1(i-1,i-1), A1(i,i-1), A1(i,i)) ;
        if(len1>0) CORE_zlarfx2(PlasmaLeft, len1 , *V1(i), conj(*TAU1(i)), A1((i-1), J1), LDX, A1(i, J1), LDX);
     }
  }
  /* **********************************************************************************************
   *   APPLY RIGHT ON THE REMAINING ELEMENT OF KERNEL 1
   * ***********************************************************************************************/
  for (i = ed; i >= st+1 ; i--){
     J1    = i+1;
     J2    = ed;
     len   = J2-J1+1;
     if(len>0){
        if(i>NB)
           /* both column (i-1) and i are on tile T2 */
           CORE_zlarfx2(PlasmaRight, len, *V1(i), *TAU1(i), A2(J1, i-1), LDX, A2(J1, i), LDX);
        else if(i==NB)
           /* column (i-1) is on tile T1 while column i is on tile T2 */
           CORE_zlarfx2(PlasmaRight, len, *V1(i), *TAU1(i), A1(J1, i-1), LDX, A2(J1, i), LDX);
        else
           /* both column (i-1) and i are on tile T1 */
           CORE_zlarfx2(PlasmaRight, len, *V1(i), *TAU1(i), A1(J1, i-1), LDX, A1(J1, i), LDX);
     }
  }
}
///////////////////////////////////////////////////////////
#undef A1
#undef A2
#undef V1
#undef TAU1
#undef V2
#undef TAU2

///////////////////////////////////////////////////////////
//                  TYPE 1-BAND Householder
///////////////////////////////////////////////////////////
//// add -1 because of C
#define A(m,n)   &(A[((m)-(n)) + LDA*((n)-1)])
#define V(m)     &(V[m-1])
#define TAU(m)   &(TAU[m-1])
static void TRD_type1bHL(int N, int NB, parsec_complex64_t *A, int LDA, parsec_complex64_t *V, parsec_complex64_t *TAU, int st, int ed) {
  int    J1, J2, len, LDX;
  int    i, IONE, ITWO;
  IONE=1;
  ITWO=2;
  (void)NB;

  if(ed <= st){
    printf("TRD_type 1bH: ERROR st and ed %d  %d \n",st,ed);
    exit(-10);
  }

  LDX = LDA-1;
  for (i = ed; i >= st+1 ; i--){
     // generate Householder to annihilate a(i+k-1,i) within the band
     *V(i)          = *A(i, (st-1));
     *A(i, (st-1))  = 0.0;
     DLARFG( &ITWO, A((i-1),(st-1)), V(i), &IONE, TAU(i) );

     // apply reflector from the left (horizontal row) and from the right for only the diagonal 2x2.
     J1  = st;
     J2  = i;
     len = J2-J1+1;
     printf("voici J1 %d J2 %d len %d \n",J1,J2,len);
     DLARFX_C('B', len , *V(i), *TAU(i), A((i-1),J1   ), LDX);
  }

  for (i = ed; i >= st+1 ; i--){
     len = min(ed,N)-i;
     if(len>0)DLARFX_C('R', len, *V(i), *TAU(i), A((i+1),(i-1)), LDX);
  }

}
#undef A
#undef V
#undef TAU
///////////////////////////////////////////////////////////

///////////////////////////////////////////////////////////
//                  TYPE 2-BAND Householder
///////////////////////////////////////////////////////////
//// add -1 because of C
#define A(m,n)   &(A[((m)-(n)) + LDA*((n)-1)])
#define V(m)     &(V[m-1])
#define TAU(m)   &(TAU[m-1])
static void TRD_type2bHL(int N, int NB, parsec_complex64_t *A, int LDA, parsec_complex64_t *V, parsec_complex64_t *TAU, int st, int ed) {
  int    J1, J2, J3, KDM2, len, LDX;
  int    i, IONE, ITWO;
  IONE=1;
  ITWO=2;

  if(ed <= st){
    printf("TRD_type 2H: ERROR st and ed %d  %d \n",st,ed);
    exit(-10);;
  }

  LDX  = LDA -1;
  KDM2 = NB-2;
  for (i = ed; i >= st+1 ; i--){
     // apply Householder from the right. and create newnnz outside the band if J3 < N
     J1  = ed+1;
     J2  = min((i+1+KDM2), N);
     J3  = min((J2+1), N);
     len = J2-J1+1;
     DLARFX_C('R', len, *V(i), *TAU(i), A(J1,(i-1)), LDX);

     // if nonzero element a(j+kd,j-1) has been created outside the band (if index < N) then eliminate it.
     len    = J3-J2; // soit 1 soit 0
     if(len>0){
         //  new nnz at TEMP=V(J3)
         *V(J3)     =            - *A(J3,(i)) * (*TAU(i)) * (*V(i));
         *A(J3,(i)) = *A(J3,(i)) + *V(J3)  * (*V(i)); //ATTENTION THIS replacement IS VALID IN FLOAT CASE NOT IN COMPLEX
         // generate Householder to annihilate a(j+kd,j-1) within the band
         DLARFG( &ITWO, A(J2,(i-1)), V(J3), &IONE, TAU(J3) );
     }
  }
               //if(id==1) return;

  for (i = ed; i >= st+1 ; i--){
     J2  = min((i+1+KDM2), N);
     J3  = min((J2+1), N);
     len    = J3-J2;
     if(len>0){
        len = min(ed,N)-i+1;
        DLARFX_C('L', len , *V(J3), *TAU(J3), A(J2, i), LDX);
     }
  }

}
#undef A
#undef V
#undef TAU
///////////////////////////////////////////////////////////

///////////////////////////////////////////////////////////
//                  TYPE 3-BAND Householder
///////////////////////////////////////////////////////////
//// add -1 because of C
#define A(m,n)   &(A[((m)-(n)) + LDA*((n)-1)])
#define V(m)     &(V[m-1])
#define TAU(m)   &(TAU[m-1])
static void TRD_type3bHL(int N, int NB, parsec_complex64_t *A, int LDA, parsec_complex64_t *V, parsec_complex64_t *TAU, int st, int ed) {
  int    J1, J2, len, LDX;
  int    i;
  (void)NB;

  if(ed <= st){
    printf("TRD_type 3H: ERROR st and ed %d  %d \n",st,ed);
    exit(-10);
  }

  LDX = LDA-1;
  for (i = ed; i >= st+1 ; i--){
     // apply rotation from the left. faire le DROT horizontal
     J1  = st;
     J2  = i;
     len = J2-J1+1;
     //len = NB+1;
     DLARFX_C('B', len , *V(i), *TAU(i), A((i-1),J1   ), LDX);
  }

  for (i = ed; i >= st+1 ; i--){
     len = min(ed,N)-i;
     if(len>0)DLARFX_C('R', len, *V(i), *TAU(i), A((i+1),(i-1)), LDX);
  }
}
#undef A
#undef V
#undef TAU
///////////////////////////////////////////////////////////

///////////////////////////////////////////////////////////
//                  grouping sched wrapper call
///////////////////////////////////////////////////////////
#if 0
int TRD_seqgralgtype(int N, int NB, parsec_complex64_t *A, int LDA, parsec_complex64_t *C, parsec_complex64_t *S, int i, int j, int m, int grsiz, int BAND) {
  int    k,shift=3;
  int    myid,colpt,stind,edind,blklastind,stepercol;
  (void) BAND;

  k   = shift/grsiz;
  stepercol =  k*grsiz == shift ? k:k+1;

   for (k = 1; k <=grsiz; k++){
      myid = (i-j)*(stepercol*grsiz) +(m-1)*grsiz + k;
      if(myid%2 ==0){
           colpt      = (myid/2)*NB+1+j-1;
           stind      = colpt-NB+1;
           edind      = min(colpt,N);
           blklastind = colpt;
           if(stind>=edind){
               printf("TRD_seqalg ERROR---------> st>=ed  %d  %d \n\n",stind, edind);
               return -10;
           }
      }else{
           colpt      = ((myid+1)/2)*NB + 1 +j -1 ;
           stind      = colpt-NB+1;
           edind      = min(colpt,N);
           if( (stind>=edind-1) && (edind==N) )
               blklastind=N;
           else
               blklastind=0;
           if(stind>=edind){
               printf("TRD_seqalg ERROR---------> st>=ed  %d  %d \n\n",stind, edind);
               return -10;

           }
      }

          if(myid == 1)
              TRD_type1bHL(N, NB, A, LDA, C, S, stind, edind);
          else if(myid%2 == 0)
              TRD_type2bHL(N, NB, A, LDA, C, S, stind, edind);
          else if(myid%2 == 1)
              TRD_type3bHL(N, NB, A, LDA, C, S, stind, edind);
          else{
              printf("COUCOU ERROR myid/2 %d\n",myid);
               return -10;
          }

      if(blklastind >= (N-1))  break;
      //if(myid==2)return;
  }   // END for k=1:grsiz
return 0;
}
///////////////////////////////////////////////////////////
#endif

///////////////////////////////////////////////////////////////////////////////////////////////////////////////
//                  REDUCTION BAND TO TRIDIAG
///////////////////////////////////////////////////////////////////////////////////////////////////////////////
#if 0
static void band_to_trd_v8seq(int N, int NB, parsec_complex64_t *A, int LDA, int INgrsiz, int INthgrsiz) {
  int myid, grsiz, shift, stt, st, ed, stind, edind, BAND;
  int blklastind, colpt;
  int stepercol,mylastid;
  parsec_complex64_t *C, *S;
  int i,j,m;
  int thgrsiz, thgrnb, thgrid, thed;
  int INFO;
  INFO=-1;
  BAND = 0;
  C   = malloc(N*sizeof(parsec_complex64_t));
  S   = malloc(N*sizeof(parsec_complex64_t));
  memset(C,0,N*sizeof(parsec_complex64_t));
  memset(S,0,N*sizeof(parsec_complex64_t));

  grsiz   = INgrsiz;
  thgrsiz = INthgrsiz;

  shift   = 3;
  if(grsiz==0)grsiz   = 6;
  if(thgrsiz==0)thgrsiz = N;

  if(LDA != (NB+1))
  {
      printf(" ERROR LDA not equal NB+1 and this code is special for LDA=NB+1. LDA=%d NB+1=%d \n",LDA,NB+1);
      return;
  }
  printf("  Version -8seq- grsiz %4d     thgrsiz %4d       N %5d      NB %5d    BAND %5d\n",grsiz,thgrsiz, N, NB, BAND);

  i = shift/grsiz;
  stepercol =  i*grsiz == shift ? i:i+1;

  i       = (N-2)/thgrsiz;
  thgrnb  = i*thgrsiz == (N-2) ? i:i+1;

  for (thgrid = 1; thgrid<=thgrnb; thgrid++){
     stt  = (thgrid-1)*thgrsiz+1;
     thed = min( (stt + thgrsiz -1), (N-2));
     for (i = stt; i <= N-2; i++){
        ed=min(i,thed);
        if(stt>ed)break;
        for (m = 1; m <=stepercol; m++){
            st=stt;
            for (j = st; j <=ed; j++){
                 myid     = (i-j)*(stepercol*grsiz) +(m-1)*grsiz + 1;
                 mylastid = myid+grsiz-1;
                 INFO = TRD_seqgralgtype(N, NB, A, LDA, C, S, i, j, m, grsiz, BAND);
                 if(INFO!=0){
                         printf("ERROR band_to_trd_v8seq INFO=%d\n",INFO);
                         return;
                 }
                 if(mylastid%2 ==0){
                     blklastind      = (mylastid/2)*NB+1+j-1;
                 }else{
                      colpt      = ((mylastid+1)/2)*NB + 1 +j -1 ;
                      stind      = colpt-NB+1;
                      edind      = min(colpt,N);
                      if( (stind>=edind-1) && (edind==N) )
                          blklastind=N;
                      else
                          blklastind=0;
                 }
                 if(blklastind >= (N-1))  stt=stt+1;
           } // END for j=st:ed
        } // END for m=1:stepercol
     } // END for i=1:N-2
  } // END for thgrid=1:thgrnb

} // END FUNCTION
///////////////////////////////////////////////////////////////////////////////////////////////////////////////
#endif

///////////////////////////////////////////////////////////////////////////////////////////////////////////////
int blgchase_ztrdv1(int NT, int N, int NB, parsec_complex64_t *A, parsec_complex64_t *V, parsec_complex64_t *TAU, int sweep, int id, int blktile) {
  int    /*edloc,*/ stloc, st, ed, KDM1, LDA;

  (void)NT;
  KDM1   = NB-1;
  LDA    = NB+1;
  /* generate the indiceslocal and global*/
  stloc  = (sweep+1)%NB;
  if(stloc==0) stloc=NB;
  /*if(id==NT-1)
     edloc = NB-1;
  else
     edloc = stloc + KDM1;
  */

  st = min(id*NB+stloc, N-1);
  ed = min(st+KDM1, N-1);
  /*
   * i     = (N-1)%ed;
   * edloc = i%NB+NB;
   */

  /* quick return in case of last tile */
  if(st==ed)
     return 0;

  /* because the kernel have been writted for fortran, add 1 */
   st = st +1;
   ed = ed +1;
  /* code for all tiles */
  if(id==blktile){
     TRD_type1bHL(N, NB, A, LDA, V, TAU, st, ed);
     TRD_type2bHL(N, NB, A, LDA, V, TAU, st, ed);
  }else{
     TRD_type3bHL(N, NB, A, LDA, V, TAU, st, ed);
     //if(id==6) return;
     TRD_type2bHL(N, NB, A, LDA, V, TAU, st, ed);
  }
  return 0;
}

#if 0
///////////////////////////////////////////////////////////////////////////////////////////////////////////////
//
///////////////////////////////////////////////////////////////////////////////////////////////////////////////
static void band_to_trd_vmpi1(int N, int NB, parsec_complex64_t *A, int LDA) {
  int NT;
  parsec_complex64_t *V, *TAU;
  int blktile, S, id, sweep;
  V   = malloc(N*sizeof(parsec_complex64_t));
  TAU = malloc(N*sizeof(parsec_complex64_t));
  memset(V,0,N*sizeof(parsec_complex64_t));
  memset(TAU,0,N*sizeof(parsec_complex64_t));

  NT = N/NB;
  if(NT*NB != N){
      printf("ERROR NT*NB not equal N \n");
      return;
  }

  //printf("voici NT, N NB %d %d %d\n",NT,N,NB);
  for (blktile = 0; blktile<NT; blktile++){
     for (S = 0; S<NB; S++){
        sweep = blktile*NB + S ;
        for (id = blktile; id<NT; id++){
               blgchase_ztrdv1 (NT, N, NB, A, V, TAU, sweep , id, blktile);
        }
     }
  }

}
///////////////////////////////////////////////////////////////////////////////////////////////////////////////
#endif

///////////////////////////////////////////////////////////////////////////////////////////////////////////////
int blgchase_ztrdv2(int NT, int N, int NB, parsec_complex64_t *A1, parsec_complex64_t *A2, parsec_complex64_t *V1, parsec_complex64_t *TAU1, parsec_complex64_t *V2, parsec_complex64_t *TAU2, int sweep, int id, int blktile) {
  int    edloc, stloc, st, ed, KDM1, LDA;

  KDM1   = NB-1;
  LDA    = NB+1;
  /* generate the indiceslocal and global*/
  stloc  = (sweep+1)%NB;
  if(stloc==0) stloc=NB;
  if(id==NT-1)
     edloc = NB-1;
  else
     edloc = stloc + KDM1;

  st = min(id*NB+stloc, N-1);
  ed = min(st+KDM1, N-1);
  /*
   * i     = (N-1)%ed;
   * edloc = i%NB+NB;
   */

  /* quick return in case of last tile */
  if(st==ed)
     return 0;

  /* code for all tiles */
  if(id==blktile){
     CORE_zhbtelr(N, NB, A1, LDA, A2, LDA, V1, TAU1, stloc, edloc);
     CORE_zhbtrce(N, NB, A1, LDA, A2, LDA, V1, TAU1, V2, TAU2, stloc, edloc, ed);
  }else{
     CORE_zhbtlrx(N, NB, A1, LDA, A2, LDA, V1, TAU1, stloc, edloc);
     CORE_zhbtrce(N, NB, A1, LDA, A2, LDA, V1, TAU1, V2, TAU2, stloc, edloc, ed);
  }
  return 0;
}

#if 0
///////////////////////////////////////////////////////////////////////////////////////////////////////////////
//
///////////////////////////////////////////////////////////////////////////////////////////////////////////////
static void band_to_trd_vmpi2(int N, int NB, parsec_complex64_t *A, int LDA) {
    int NT;
    parsec_complex64_t *V, *TAU;
    int blktile, S, id, sweep;
    V   = malloc(N*sizeof(parsec_complex64_t));
    TAU = malloc(N*sizeof(parsec_complex64_t));
    memset(V,0,N*sizeof(parsec_complex64_t));
    memset(TAU,0,N*sizeof(parsec_complex64_t));

    NT = N/NB;
    if(NT*NB != N){
        printf("ERROR NT*NB not equal N \n");
        return;
    }

    //printf("voici NT, N NB %d %d %d\n",NT,N,NB);
    for (blktile = 0; blktile<NT; blktile++){
        for (S = 0; S<NB; S++){
            sweep = blktile*NB + S ;
            for (id = blktile; id<NT; id++){
                //printf("voici  blktile %d    S %d     id %d   sweep %d   \n",blktile, S, id, sweep);
                blgchase_ztrdv2 (NT, N, NB,
                                A+(id*NB*LDA), A+((id+1)*NB*LDA),
                                V+(id*NB), TAU+(id*NB),
                                V+((id+1)*NB), TAU+((id+1)*NB),
                                sweep, id, blktile);
            }
            //if(sweep==0) return;
        }
    }

}
///////////////////////////////////////////////////////////////////////////////////////////////////////////////
#endif


text:
/*
** Various utilities for LISA algorithm
**
** G.Lohmann, Jan 2017
*/
#include <viaio/VImage.h>
#include <viaio/Vlib.h>
#include <viaio/mu.h>
#include <viaio/option.h>

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

#include <gsl/gsl_math.h>
#include <gsl/gsl_histogram.h>

#define SQR(x) ((x)*(x))
#define ABS(x) ((x) > 0 ? (x) : -(x))

extern float kth_smallest(float *a, size_t n, size_t k);
#define Median(a,n) kth_smallest(a,n,(((n)&1)?((n)/2):(((n)/2)-1)))

/* update histogram */
void HistoUpdate(VImage src1,gsl_histogram *hist)
{
  float u,tiny = 1.0e-6;
  size_t i;
  float xmin = gsl_histogram_min (hist);
  float xmax = gsl_histogram_max (hist);

  float *pp1 = VImageData(src1);
  for (i=0; i<VImageNPixels(src1); i++) {
    u = *pp1++;
    if (ABS(u) < tiny) continue;
    if (u > xmax) u = xmax-tiny;
    if (u < xmin) u = xmin+tiny;
    gsl_histogram_increment (hist,u);
  }
}

/* In moderately skewed or asymmetrical distribution (Pearson) */
float VGetMode(VImage src)
{
  size_t i,n=0;
  float u,tiny=1.0e-8;
  VFloat *pp = VImageData(src);
  for (i=0; i<VImageNPixels(src); i++) {
    u = (*pp++);
    if (fabs(u) > tiny) n++;
  }
  float *data = (float *) VCalloc(n,sizeof(float));

  double sum=0,nx=0;
  n=0;
  pp = VImageData(src);
  for (i=0; i<VImageNPixels(src); i++) {
    u = (*pp++);
    if (fabs(u) > tiny) {
      data[n] = u;
      sum += u;
      n++;
    }
  }
  nx = (double)n;
  float median = Median(data,n);
  float mean = (float)(sum/nx);
  float mode = 3.0*median - 2.0*mean;
  /* fprintf(stderr," mean: %f,  median: %f,  mode: %f\n",mean,median,mode); */
  VFree(data);
  return mode;
}

/* scale z-values */
void VZScale(VImage src,float mode,float stddev)
{
  size_t i=0;
  float u=0,tiny=1.0e-8;
  VFloat *pp=VImageData(src);
  for (i=0; i<VImageNPixels(src); i++) {
    u = (*pp);
    if (fabs(u) > tiny) (*pp) = (u-mode)/stddev;
    pp++;
  }
}

/* get max z-value */
float VZMax(VImage src)
{
  size_t i=0;
  float u=0,umax=0;
  VFloat *pp=VImageData(src);
  for (i=0; i<VImageNPixels(src); i++) {
    u = (*pp++);
    if (u > umax) umax = u;
  }
  return umax;
}

/* get image variance */
double VImageVar(VImage src)
{
  size_t i=0;
  double u=0,s1=0,s2=0,nx=0,tiny=1.0e-8;
  VFloat *pp=VImageData(src);
  for (i=0; i<VImageNPixels(src); i++) {
    u = (double)(*pp++);
    if (ABS(u) < tiny) continue;
    s1 += u;
    s2 += u*u;
    nx++;
  }
  if (nx < 3.0) VError(" number of nonzero voxels: %g\n",nx);
  double mean = s1/nx;
  double var = (s2 - nx * mean * mean) / (nx - 1.0);
  return var;
}

/* count nonzero voxels */
void VImageCount(VImage src)
{
  float u=0;
  size_t i=0,npos=0;
  VFloat *pp=VImageData(src);
  for (i=0; i<VImageNPixels(src); i++) {
    u = (float)(*pp++);
    if (u > 0) npos++;
  }
}

/* check histogram range, image is normalized so [-10,10] should be okay */
void VGetHistRange(VImage src,double *hmin,double *hmax)
{
  size_t i=0;
  double u=0,x=0,y=0;
  double zmin = -10.0;
  double zmax = 10.0;

  VFloat *pp=VImageData(src);
  size_t npix=0,npos=0,nneg=0;
  for (i=0; i<VImageNPixels(src); i++) {
    u = (*pp++);
    if (u < zmin) nneg++;
    if (u > zmax) npos++;
    if (fabs(u) > 0) npix++;
  }
  if (npix < 1) VError(" image has no non-zero pixels");

  if (npos > 0 || nneg > 0) {
    x = (double)npos/(double)npix;
    y = (double)nneg/(double)npix;
    if (x > 0.01 || y > 0.01) {
      zmin = -20.0;
      zmax = 20.0;
    }
  }

  /* output */
  *hmin = zmin;
  *hmax = zmax;
}

/* remove isolated voxels */
void VIsolatedVoxels(VImage src,float threshold)
{
  int b,r,c,bb,rr,cc;
  float u=0;
  size_t i;

  VImage tmp = VCreateImage(VImageNBands(src),VImageNRows(src),VImageNColumns(src),VBitRepn);
  VFillImage(tmp,VAllBands,0);

  VFloat *pp=VImageData(src);
  VBit *pa=VImageData(tmp);
  for (i=0; i<VImageNPixels(src); i++) {
    u = (*pp++);
    if (u > threshold) *pa = 1;
    pa++;
  }

  for (b=1; b<VImageNBands(src)-1; b++) {
    for (r=1; r<VImageNRows(src)-1; r++) {
      for (c=1; c<VImageNColumns(src)-1; c++) {

	int n=0;
	for (bb=b-1; bb<=b+1; bb++) {
	  for (rr=r-1; rr<=r+1; rr++) {
	    for (cc=c-1; cc<=c+1; cc++) {
	      if (VPixel(tmp,bb,rr,cc,VBit) > 0) n++;
	      if (n >= 2) goto skip;
	    }
	  }
	}
	if (n < 2) VPixel(src,b,r,c,VFloat) = 0;
      skip: ;
      }
    }
  }
  VDestroyImage(tmp);
}

ext: c

max_stars_repo_name: zrajna/lipsia

max_issues_repo_name: zrajna/lipsia

max_forks_repo_name: zrajna/lipsia


text:
#include <gsl/gsl_sf_airy.h>
#include <gsl/gsl_errno.h>

/* function names generated by removing  "gsl_sf" from the beginning
   of the name.  Thus gsl_sf_airy_Ai goes to airy_Ai. */

static int sf_mode[] = { GSL_PREC_DOUBLE, GSL_PREC_SINGLE, GSL_PREC_APPROX } ;

void airy_Ai_e(double *x, int *len, int *mode, double *val, double *err, int *status)
{
	int i;
	gsl_sf_result result;
	gsl_set_error_handler_off();

	for(i = 0; i< *len ; i++){
		status[i] = gsl_sf_airy_Ai_e(x[i], sf_mode[*mode], &result) ;
		val[i] = result.val;
		err[i] = result.err;
	}
}

void airy_Bi_e(double *x, int *len, int *mode, double *val, double *err, int *status)
{
	int i;
	gsl_sf_result result;

	for(i = 0; i< *len ; i++){
		status[i] = gsl_sf_airy_Bi_e(x[i], sf_mode[*mode], &result) ;
		val[i] = result.val;
		err[i] = result.err;
	}
}

void airy_Ai_scaled_e(double *x, int *len, int *mode, double *val, double *err, int *status)
{
	int i;
	gsl_sf_result result;

	for(i = 0; i< *len ; i++){
		status[i] = gsl_sf_airy_Ai_scaled_e(x[i], sf_mode[*mode], &result) ;
		val[i] = result.val;
		err[i] = result.err;
	}
}

void airy_Bi_scaled_e(double *x, int *len, int *mode, double *val, double *err, int *status)
{
	int i;
	gsl_sf_result result;

	for(i = 0; i< *len ; i++){
		status[i] = gsl_sf_airy_Bi_scaled_e(x[i], sf_mode[*mode], &result) ;
		val[i] = result.val;
		err[i] = result.err;
	}
}

void airy_Ai_deriv_e(double *x, int *len, int *mode, double *val, double *err, int *status)
{
	int i;
	gsl_sf_result result;

	for(i = 0; i< *len ; i++){
		status[i] = gsl_sf_airy_Ai_deriv_e(x[i], sf_mode[*mode], &result) ;
		val[i] = result.val;
		err[i] = result.err;
	}
}

void airy_Bi_deriv_e(double *x, int *len, int *mode, double *val, double *err, int *status)
{
	int i;
	gsl_sf_result result;

	for(i = 0; i< *len ; i++){
		status[i] = gsl_sf_airy_Bi_deriv_e(x[i], sf_mode[*mode], &result) ;
		val[i] = result.val;
		err[i] = result.err;
	}
}

void airy_Ai_deriv_scaled_e(double *x, int *len, int *mode, double *val, double *err, int *status)
{
	int i;
	gsl_sf_result result;

	for(i = 0; i< *len ; i++){
		status[i] = gsl_sf_airy_Ai_deriv_scaled_e(x[i], sf_mode[*mode], &result) ;
		val[i] = result.val;
		err[i] = result.err;
	}
}

void airy_Bi_deriv_scaled_e(double *x, int *len, int *mode, double *val, double *err, int *status)
{
	int i;
	gsl_sf_result result;

	for(i = 0; i< *len ; i++){
		status[i] = gsl_sf_airy_Bi_deriv_scaled_e(x[i], sf_mode[*mode], &result) ;
		val[i] = result.val;
		err[i] = result.err;
	}
}

void airy_zero_Ai_e(int *n, int *len, double *val, double *err, int *status)
{
	int i;
	gsl_sf_result result;
	for(i = 0; i< *len ; i++){
		if(n[i]>0){
			status[i] = gsl_sf_airy_zero_Ai_e(n[i], &result) ;
		} else {
			result.val=0.0;
			result.err=GSL_EDOM;
		}
		val[i] = result.val;
		err[i] = result.err;
	}
}

void airy_zero_Bi_e(int *n, int *len, double *val, double *err, int *status)
{
	int i;
	gsl_sf_result result;
	for(i = 0; i< *len ; i++){
		if(n[i]>0){
			status[i] = gsl_sf_airy_zero_Bi_e(n[i], &result) ;
		} else {
			result.val=0.0;
			result.err=GSL_EDOM;
		}
		val[i] = result.val;
		err[i] = result.err;
	}
}

void airy_zero_Ai_deriv_e(int *n, int *len, double *val, double *err, int *status)
{
	int i;
	gsl_sf_result result;
	for(i = 0; i< *len ; i++){
		if(n[i]>0){
			status[i] = gsl_sf_airy_zero_Ai_deriv_e(n[i], &result) ;
		} else {
			result.val=0.0;
			result.err=GSL_EDOM;
		}
		val[i] = result.val;
		err[i] = result.err;
	}
}

void airy_zero_Bi_deriv_e(int *n, int *len, double *val, double *err, int *status)
{
	int i;
	gsl_sf_result result;
	for(i = 0; i< *len ; i++){
		if(n[i]>0){
			status[i] = gsl_sf_airy_zero_Bi_deriv_e(n[i], &result) ;
		} else {
			result.val=0.0;
			result.err=GSL_EDOM;
		}
		val[i] = result.val;
		err[i] = result.err;
	}
}

max_stars_repo_name: kbroman/fluctuationDomains

max_issues_repo_name: kbroman/fluctuationDomains

max_forks_repo_name: kbroman/fluctuationDomains


text:
#ifndef __calcAinvb__
#define __calcAinvb__

#include <gsl/gsl_linalg.h>
#include <gsl/gsl_matrix_double.h>
//#include <gsl/gsl_multimin.h>
#include <gsl/gsl_multiroots.h>
//#include "my_misc.c"
//#define GSL 1
#define GSL 1
int GSLtag;
void least_gsl_error_handler( const char * reason,
			      const char * file,
			      int line,
			      int gsl_errno)
{
  //  printf("%s:%d: %s (gsl_error: %s)\n",
  //	 file, line, reason, gsl_strerror(gsl_errno));
  printf("%s:%d: %s (gsl_error %d)\n",
	 file, line, reason,gsl_errno);
  GSLtag=1;
}
double VarY;
void calc_Ainvb(double *M, double *A_data, double *b_data, int nx, int ndata,int nxmarume)
{
  /* b= A M  --> M= Ainv b
     b  in R^{ndata x 1}
     A  in R^{ndata x nx}   (ndata>nx)
     M  in R^{nx    x 1}
   */
  int i;
  gsl_vector *S;
  gsl_vector *work;
  gsl_vector *x;
  gsl_vector *b;
  gsl_matrix *V;
  gsl_matrix *A;
  if(ndata<=0){
    M[0]=1;
    for(i=1;i<nx;i++) M[i]=1./nx;
  }
  else{
    A =gsl_matrix_alloc(ndata,nx);
    b =gsl_vector_alloc(ndata);b->data=b_data;
    A->data=A_data;
    //  for(i=0;i<ndata;i++){printf("%e %e %d #bi,Ai\n",A->data[i],b->data[i],i); }

    V = gsl_matrix_alloc (nx,nx);
    S = gsl_vector_alloc (nx);
    work = gsl_vector_alloc(nx);
    x = gsl_vector_alloc(nx);

    gsl_set_error_handler( &least_gsl_error_handler );GSLtag=0;

#if GSL == 1
    gsl_linalg_SV_decomp_jacobi (A, V, S); //higer accuracy than Golub-Reinsh
#elif GSL == 2
    gsl_linalg_SV_decomp (A, V, S, work); //Golub-Reinsh
#elif GSL == 3
    {//faster
      gsl_matrix *X;
      X = gsl_matrix_alloc (nx,nx);
      gsl_linalg_SV_decomp_mod (A, X, V, S, work);
      gsl_matrix_free(X);
    }
#endif
    //    for(i=nxmarume;i<nx;i++){

    fprintf(stdout,"Lambda");
    for(i=0;i<nx;i++){
      //      if(nx>40) break;
      //      fprintf(stdout,"S->data[%d]=%e;",i,S->data[i]);
      fprintf(stdout,"(%d)%e",i,S->data[i]);
    }
    fprintf(stdout,"\n");
    VarY=0;
    if(GSLtag==1){
      fprintf(stderr,"\nEigenValue=");for(i=0;i<nx;i++) fprintf(stderr,"%e ",S->data[i]); fprintf(stderr,"\n");
      for(i=0;i<nx;i++) M[i]=1./nx;
    }
    else{
      if(nxmarume>=0){//Principal Component Analysis
	for(i=nxmarume;i<nx;i++) S->data[i]=0;
	for(i=0;i<nxmarume;i++) VarY+=S->data[i];
	gsl_linalg_SV_solve(A,V,S,b,x);
	for(i=0;i<nx;i++) M[i]=gsl_vector_get(x,i);
	//    {//check
	//      double err=0,yhat;
	//      int j;
	//      for(i=0;i<ndata;i++){
	//	if(i>10) break;
	//	yhat=b_data[i];
	//	for(j=0;j<nx;j++){
	//	  yhat-= M[j]*A_data[i*nx+j];
	//	}
	//	err+=square(yhat);
	//	fprintf(stderr,"err=%e=%e/%d. yhat=%e=(%e",err/(i+1.),err,i,yhat,b_data[i]);
	//	for(j=0;j<nx;j++) fprintf(stderr,"-(%e)*(%e)",M[j],A_data[i*nx+j]);
	//	fprintf(stderr,")^2\n");
	//      }
	//    }
      }
      else{//Minor Component Analysis
	int ii=0;
	int im=-nxmarume;
	for(i=nx-1;i>=0;i--){
	  if(S->data[i]>1e-20){
	    ii++;
	    VarY+=S->data[i];
	    if(ii>=im) break;
	  }
	}
	int j;

	double M0=-V->data[i];
	for(j=1;j<nx;j++){
	  M[j-1]=V->data[j*nx+i]/M0;
	}
//	double M0=-V->data[i*nx];
//	for(j=1;j<nx;j++){
//	  M[j-1]=V->data[i*nx+j]/M0;
//	}

      }
    }
    gsl_vector_free(S);
    gsl_vector_free(work);
    gsl_vector_free(x);
    gsl_matrix_free(V);
    gsl_matrix_free(A);
    gsl_vector_free(b);
  }
}
void calc_AtinvbMCA(double *M, double *At_data, double *b_data, int nx, int ndata,int nxmarume)
{
  /* b= A M  --> M= Ainv b
     b  in R^{ndata x 1}
     A=At  in R^{ndata x nx}
     M  in R^{nx    x 1}
   */
  int nx1=nx+1,i,j,jnx1;
  gsl_matrix *A  =gsl_matrix_alloc(ndata,nx1);
  for(j=0;j<ndata;j++){
    jnx1=j*nx1;
    A->data[jnx1]=b_data[j];
    for(i=0;i<nx;i++){
      A->data[jnx1+i+1]=At_data[i*ndata+j];
    }
  }
  calc_Ainvb(M, (double *)(A->data), b_data, nx1, ndata, -nxmarume);

  gsl_matrix_free(A);
}
void calc_Atinvb(double *M, double *At_data, double *b_data, int nx, int ndata,int nxmarume)
{
  /* b= A M  --> M= Ainv b
     b  in R^{ndata x 1}
     A=At  in R^{ndata x nx}
     M  in R^{nx    x 1}
   */
  gsl_matrix *At =gsl_matrix_alloc(nx,ndata);
  gsl_matrix *A  =gsl_matrix_alloc(ndata,nx);
  At->data=At_data;
  gsl_matrix_transpose_memcpy(A, At);
  calc_Ainvb(M, (double *)(A->data), b_data, nx, ndata, nxmarume);
  gsl_matrix_free(A);
  gsl_matrix_free(At);
}

#endif // __calcAinvb__


text:
/** @file */

#ifndef __CCL_F3D_H_INCLUDED__
#define __CCL_F3D_H_INCLUDED__

#include <gsl/gsl_spline.h>
#include <gsl/gsl_interp2d.h>
#include <gsl/gsl_spline2d.h>

CCL_BEGIN_DECLS

/**
 * Struct for accelerated linear interpolation.
 */
typedef struct {
  int ia_last; /**< Last index found */
  double amin; /**< Minimum a-value within range */
  double amax; /**< Maximum a-value within range */
  int na; /**< Number of a-values held */
  double *a_arr; /**< Array of a-values */
} ccl_a_finder;

/**
 * Creates a new ccl_a_finder structure from an array
 * of scale factors.
 * @param na Number of elements held by a_arr
 * @param a_arr array of scale factors over which linear interpolation will be carried out.
 */
ccl_a_finder *ccl_a_finder_new(int na, double *a_arr);

/**
 * ccl_a_finder destructor.
 */
void ccl_a_finder_free(ccl_a_finder *finda);

/**
 * Find index corresponding to scale factor value a
 * such that finda->a_arr[index]<a<finda->a_arr[index+1].
 * @param finda ccl_a_finder.
 * @param a scale factor value.
 */
int ccl_find_a_index(ccl_a_finder *finda, double a);

/**
 * Struct containing a 3D trispectrum
 */
typedef struct {
  double lkmin,lkmax; /**< Edges in log(k)*/
  int na; /**< Number of a values */
  double *a_arr; /**< Array of a values at which this is sampled */
  int is_product; /**< Is this factorizable as f(k1,a)*g(k2,a)? */
  int extrap_order_lok; /**< Order of extrapolating polynomial in log(k) for low k (0, 1 or 2)*/
  int extrap_order_hik; /**< Order of extrapolating polynomial in log(k) for high k (0, 1 or 2)*/
  ccl_f2d_extrap_growth_t extrap_linear_growth;  /**< Extrapolation type at high redshifts*/
  int is_log; /**< Do I hold the values of log(f(k,a))?*/
  double growth_factor_0; /**< Constant extrapolating growth factor*/
  int growth_exponent; /**< Power to which growth should be exponentiated*/
  ccl_f2d_t *fka_1; /**< If is_product=True, then this holds the first factor f(k,a) */
  ccl_f2d_t *fka_2; /**< If is_product=True, then this holds the second factor g(k,a) */
  gsl_spline2d **tkka; /**< Array of 2D (k1,k2) splines (one for each value of a). */
} ccl_f3d_t;

/**
 * Create a ccl_f3d_t structure.
 * @param na number of elements in a_arr.
 * @param a_arr array of scale factor values at which the function is defined. The array should be ordered.
 * @param nk number of elements of lk_arr.
 * @param lk_arr array of logarithmic wavenumbers at which the function is defined (i.e. this array contains ln(k), NOT k). The array should be ordered.
 * @param tkka_arr array of size na * nk * nk containing the 3D function. The 3D ordering is such that fka_arr[ik1+nk*(ik2+nk*ia)] = f(k1=exp(lk_arr[ik1]),k2=exp(lk_arr[ik2],a=a_arr[ia]).
 * @param fka1_arr array of size nk * na containing the first factor f1 making up the total function if it's factorizable such that f(k1,k2,a) = f1(k1,a)*f2(k2,a). The 2D ordering of this array should be such that fka1_arr[ik+nk*ia] = f1(k=exp(lk_arr[ik]),a=a_arr[ia]). Only relevant if is_product is true.
 * @param fka2_arr same as fka1_arr for the second factor.
 * @param is_product if not 0, fka1_arr and fka2_arr will be used as 2-D arrays to construct a factorizable 3D function f(k1,k1,a) = f1(k1,a)*f2(k2,a).
 * @param extrap_order_lok Order of the polynomial that extrapolates on wavenumbers smaller than the minimum of lk_arr. Allowed values: 0 (constant) and 1 (linear extrapolation). Extrapolation happens in ln(k).
 * @param extrap_order_hik Order of the polynomial that extrapolates on wavenumbers larger than the maximum of lk_arr. Allowed values: 0 (constant) and 1 (linear extrapolation). Extrapolation happens in ln(k).
 * @param extrap_linear_growth: ccl_f2d_extrap_growth_t value defining how the function with scale factors below the interpolation range. Allowed values: ccl_f2d_cclgrowth (scale with the CCL linear growth factor), ccl_f2d_constantgrowth (scale by multiplying the function at the earliest available scale factor by a constant number, defined by `growth_factor_0`), ccl_f2d_no_extrapol (throw an error if the function is ever evaluated outside the interpolation range in a). Note that, above the interpolation range (i.e. for low redshifts), the function will be assumed constant.
 * @param is_tkka_log: if not zero, `tkka_arr` contains ln(f(k1,k2,a)) instead of f(k1,k2,a) (and likewise for fka1_arr and fka2_arr).
 * @param growth_factor_0: growth factor outside the range of scale factors held by a_arr. Irrelevant if extrap_linear_growth!=ccl_f2d_constantgrowth.
 * @param growth_exponent: power to which the extrapolating growth factor should be exponentiated when extrapolating (e.g. usually 4 for trispectra).
 * @param interp_type: 2D interpolation method in k1,k2 space. Currently only ccl_f2d_3 is implemented (bicubic interpolation). Note that linear interpolation is used between values of the scale factor.
 * @param status Status flag. 0 if there are no errors, nonzero otherwise.
 */
ccl_f3d_t *ccl_f3d_t_new(int na,double *a_arr,
			 int nk,double *lk_arr,
                         double *tkka_arr,
			 double *fka1_arr,
			 double *fka2_arr,
			 int is_product,
			 int extrap_order_lok,
			 int extrap_order_hik,
			 ccl_f2d_extrap_growth_t extrap_linear_growth,
			 int is_tkka_log,
			 double growth_factor_0,
			 int growth_exponent,
			 ccl_f2d_interp_t interp_type,
			 int *status);

/**
 * Evaluate 3D function of k1, k2 and a defined by ccl_f3d_t structure.
 * @param f3d ccl_f3d_t structure defining f(k1,k2,a).
 * @param lk1 Natural logarithm of the wavenumber.
 * @param lk2 Natural logarithm of the wavenumber.
 * @param a Scale factor.
 * @param finda Helper structure used to accelerate the scale factor interpolation.
 * @param cosmo ccl_cosmology structure, only needed if evaluating f(k1,k2,a) at small scale factors outside the interpolation range, and if fka was initialized with extrap_linear_growth = ccl_f2d_cclgrowth.
 * @param status Status flag. 0 if there are no errors, nonzero otherwise.
 */
double ccl_f3d_t_eval(ccl_f3d_t *f3d,double lk1,double lk2,double a,ccl_a_finder *finda,
                      void *cosmo, int *status);

/**
 * F3D structure destructor.
 * Frees up all memory associated with a f3d structure.
 * @param f3d Structure to be freed.
 */
void ccl_f3d_t_free(ccl_f3d_t *f3d);

/**
 * Make a copy of a ccl_f3d_t structure.
 * @param f3d_o old ccl_f3d_t structure.
 * @param status Status flag. 0 if there are no errors, nonzero otherwise.
 */
ccl_f3d_t *ccl_f3d_t_copy(ccl_f3d_t *f3d_o, int *status);

/**
 * Create a ccl_a_finder from the array of scale factors held
 * by a ccl_f3d_t structure.
 * @param f3d ccl_f3d_t structure.
 */
ccl_a_finder *ccl_a_finder_new_from_f3d(ccl_f3d_t *f3d);

CCL_END_DECLS

#endif

ext: h

max_stars_repo_name: Jappenn/CCL

max_issues_repo_name: Jappenn/CCL

max_forks_repo_name: Jappenn/CCL


text:
#include <errno.h>
#include <getopt.h>
#include <gsl/gsl_spmatrix.h>
#include <stdio.h>
#include <stdlib.h>
#include <sys/stat.h>

static int verbose = 0;
static int help = 0;

int test (size_t m,
	  size_t n,
	  size_t r,
	  size_t c,
	  size_t nnz,
	  const size_t *ind,
	  const size_t *ptr,
	  const double *data,
          int trials,
          int verbose);

static void usage () {
  fprintf(stderr,"usage: spmv [options] <input>\n"
  "  <input>                    MatrixMarket file (multiply this matrix)\n"
  "  -r, --block_r <arg>        Row block size\n"
  "  -c, --block_c <arg>        Column block size\n"
  "  -t, --trials <arg>         Number of trials to run\n"
  "  -v, --verbose              Verbose mode\n"
  "  -q, --quiet                Quiet mode\n"
  "  -h, --help                 Display help message\n");
}

int main (int argc, char **argv) {

  size_t b_r = 1;
  size_t b_c = 1;
  int trials = 1;

  /* Beware. Option parsing below. */
  long longarg;
  double doublearg;
  while (1) {
    static char *options = "r:c:t:vqh";
    static struct option long_options[] = {
        {"block_r",  required_argument, 0, 'r'},
        {"block_c",  required_argument, 0, 'c'},
        {"trials",   required_argument, 0, 't'},
        {"verbose",   no_argument, &verbose, 1},
        {"quiet",     no_argument, &verbose, 0},
        {"help",      no_argument, &help,    1},
        {0, 0, 0, 0}
      };

    /* getopt_long stores the option index here. */
    int option_index = 0;

    int c = getopt_long (argc, argv, options,
                     long_options, &option_index);

    /* Detect the end of the options. */
    if (c == -1)
      break;

    if (c == 0 && long_options[option_index].flag == 0)
      c = long_options[option_index].val;

    switch (c) {
      case 0:
        /* If this option set a flag, do nothing else now. */
        break;

      case 'r':
        errno = 0;
        longarg = strtol(optarg, 0, 10);
        if (errno != 0 || longarg < 1) {
          printf("option -r takes an integer column block size >= 1\n");
          usage();
          return 1;
        }
        b_r = longarg;
        break;

      case 'c':
        errno = 0;
        longarg = strtol(optarg, 0, 10);
        if (errno != 0 || longarg < 1) {
          printf("option -c takes an integer column block size >= 1\n");
          usage();
          return 1;
        }
        b_c = longarg;
        break;

      case 't':
        errno = 0;
        longarg = strtol(optarg, 0, 10);
        if (errno != 0 || longarg < 1) {
          printf("option -t takes an integer number of trials >= 1\n");
          usage();
          return 1;
        }
        trials = longarg;
        break;

      case 'v':
        verbose = 1;
        break;

      case 'q':
        verbose = 0;
        break;

      case 'h':
        help = 1;
        break;

      case '?':
        usage();
        return 1;

      default:
        abort();
    }
  }

  if (help) {
    printf("Run a fill estimation algorithm!\n");
    usage();
    return 0;
  }

  if (argc - optind > 1) {
    printf("<input> cannot be more than one file\n");
    usage();
    return 1;
  }

  if (argc - optind < 1) {
    printf("<input> not specified\n");
    usage();
    return 1;
  }

  struct stat statthing;
  if (stat(argv[optind], &statthing) < 0 || !S_ISREG(statthing.st_mode)){
    printf("<input> must be filename of MatrixMarket matrix\n");
    usage();
    return 1;
  }

  FILE *f = fopen(argv[optind], "r");
  gsl_spmatrix *triples = gsl_spmatrix_fscanf(f);
  fclose(f);
  if (triples == 0) {
    printf("<input> must be filename of MatrixMarket matrix\n");
    usage();
    return 1;
  }

  int ret = test(triples->size1, triples->size2, b_r, b_c, triples->nz, triples->i, triples->p, triples->data, trials, verbose);

  gsl_spmatrix_free(triples);

  return ret;

}


text:
#include <stdlib.h>
#include <stdio.h>
#include <stdint.h>
#include <dc1394/dc1394.h>
#include <math.h>
#include <fitsio.h>
#include <gsl/gsl_statistics.h>
#include <xpa.h>
#include <sys/time.h>
#include <string.h>
#include <inttypes.h>

#define NXPA 10

typedef struct _Box {
    double  x;
    double  y;
    double  fwhm;
    double  cenx;
    double  ceny;
    double  counts;
    double  background;
    double  noise;
    double  sigmaxythresh;
    double  sigmaxy;
    double  sigmafwhmthresh;
    double  sigmafwhm;
    int     r;
} Box;

typedef struct _Back {
    double x;
    double y;
    double background;
    double sigma;
    int r;
    int width;
} Back;

Box box[3];
Back back;
long nelements, naxes[2], fpixel;
int boxsize = 30;
double pixel_scale = 1.22;

double stardist(int i, int j) {
    return( sqrt( (box[i].cenx-box[j].cenx)*(box[i].cenx-box[j].cenx) +
		  (box[i].ceny-box[j].ceny)*(box[i].ceny-box[j].ceny) ) );
}

/* this routine uses a. tokovinin's modified equation given in
   2002, PASP, 114, 1156
*/
double seeing(double var, double d, double r) {
    double lambda;
    double b, K, seeing;

    lambda = 0.65e-6;
    b = r/d;

    /* pixel scale in "/pixel, convert to rad */
    var = var*pow(pixel_scale/206265.0,2);

    K = 0.364*(1.0 - 0.532*pow(b, -1/3) - 0.024*pow(b, -7/3));

    seeing = 206265.0*0.98*pow(d/lambda, 0.2)*pow(var/K, 0.6);

    return seeing;
}

/* this routine uses the classic DIMM equation */
double old_seeing(double var, double d, double r) {
    double lambda;
    double r0;

    lambda = 0.65e-6;

    /* pixel scale in "/pixel, convert to rad */
    var = var*pow(pixel_scale/206265.0,2);

    r0 = pow(2.0*(lambda*lambda)*(
		 ( 0.1790*pow(d, (-1.0/3.0)) -
		   0.0968*pow(r, (-1.0/3.0)) )/var
		 ), 0.6);

    return 206265.0*0.98*lambda/r0;
}

/* measure the background in an annulus around the spot pattern */
int background(char *image, int imwidth, int imheight) {

    int i, j, backpix;
    int low_y, up_y, low_x, up_x;
    double dist, sum, sumsq;

    backpix = 0;
    sum = 0.0;
    sumsq = 0.0;

    low_y = back.y - back.r - back.width;
    up_y  = back.y + back.r + back.width;
    low_x = back.x - back.r - back.width;
    up_x  = back.x + back.r + back.width;
    if (low_y < 0) {
	low_y = 0;
    }
    if (up_y >= imheight) {
	up_y = imheight;
    }
    if (low_x < 0) {
	low_x = 0;
    }
    if (up_x >= imwidth) {
	up_x = imwidth;
    }

    for (i=low_y; i<up_y; i++) {
	for (j=low_x; j<up_x; j++) {
	    dist = sqrt(pow(back.x-j, 2) + pow(back.y-i, 2));
	    if (dist >= back.r && dist <= back.r+back.width) {
		sum += image[i*naxes[0]+j];
		backpix++;
	    }
	}
    }

    back.background = sum/backpix;
    //back.background = 0.0;

    for (i=low_y; i<up_y; i++) {
	for (j=low_x; j<up_x; j++) {
	    dist = sqrt(pow(back.x-j, 2) + pow(back.y-i, 2));
	    if (dist >= back.r && dist <= back.r+back.width) {
		sumsq += (image[i*naxes[0]+j]-back.background)*
		    (image[i*naxes[0]+j]-back.background);
	    }
	}
    }

    back.sigma = sqrt(sumsq/backpix);
    //back.sigma = 1.0;

    return 1;

}

/* measure centroid using center-of-mass algorithm */
int centroid(char *image, int imwidth, int imheight, int num) {

    int i, j;
    double  sum    = 0.0;
    double  sumx   = 0.0;
    double  sumxx  = 0.0;
    double  sumy   = 0.0;
    double  sumyy  = 0.0;
    double  val = 0.0;
    double  gain = 0.5;
    double  rmom;
    double  dist;
    double nsigma = 5.0;
    int low_y, up_y, low_x, up_x;
    int sourcepix = 0;

    low_y = box[num].y - box[num].r;
    up_y  = box[num].y + box[num].r;
    low_x = box[num].x - box[num].r;
    up_x  = box[num].x + box[num].r;
    if (low_y < 0) {
	low_y = 0;
    }
    if (up_y >= imheight) {
	up_y = imheight;
    }
    if (low_x < 0) {
	low_x = 0;
    }
    if (up_x >= imwidth) {
	up_x = imwidth;
    }

    for (i=low_y; i<up_y; i++) {
	for (j=low_x; j<up_x; j++) {

	    dist = sqrt(pow(box[num].x-j, 2) + pow(box[num].y-i, 2));
	    if (dist <= box[num].r) {
		val = image[i*naxes[0]+j] - back.background;
		if (val >= nsigma*back.sigma) {
		    sum   += val;
		    sumx  += val*j;
		    sumxx += val*j*j;
		    sumy  += val*i;
		    sumyy += val*i*i;
		    sourcepix++;
		}
	    }

	}
    }

    if ( sum <= 0.0 ) {
	box[num].sigmaxy = -1.0;
	box[num].sigmafwhm = -1.0;
	//box[num].x = imwidth/2.0;
	//box[num].y = imheight/2.0;
	box[num].fwhm = -1.0;
    } else {
	rmom = ( sumxx - sumx * sumx / sum + sumyy - sumy * sumy / sum ) / sum;

	if ( rmom <= 0 ) {
	    box[num].fwhm = -1.0;
	} else {
	    box[num].fwhm = sqrt(rmom)  * 2.354 / sqrt(2.0);
	}

	box[num].counts = sum;
	box[num].cenx   = sumx / sum;
	box[num].ceny   = sumy / sum;
	box[num].x += gain*(box[num].cenx - box[num].x);
	box[num].y += gain*(box[num].ceny - box[num].y);
	box[num].sigmaxy= box[num].noise * sourcepix / box[num].counts / sqrt(6.0);
	box[num].sigmafwhm = box[num].noise * pow(sourcepix,1.5) / 10.
	    / box[num].fwhm / box[num].counts
	    * 2.354 * 2.354 / 2.0;
    }

    return 1;

}

int grab_frame(dc1394camera_t *cam, char *buf, int nbytes) {
    dc1394video_frame_t *frame=NULL;
    dc1394error_t err;

    err = dc1394_capture_dequeue(cam, DC1394_CAPTURE_POLICY_WAIT, &frame);
    if (err != DC1394_SUCCESS) {
	dc1394_log_error("Unable to capture.");
	dc1394_capture_stop(cam);
	dc1394_camera_free(cam);
	exit(1);
    }

    memcpy(buf, frame->image, nbytes);
    dc1394_capture_enqueue(cam, frame);
    return 1;
}

int add_gaussian(char *buffer, float cenx, float ceny, float a, float sigma) {
    float gauss, rsq;
    int i, j, low_x, up_x, low_y, up_y, size;

    size = 30;
    low_x = (int)(cenx-size);
    up_x = (int)(cenx+size);
    low_y = (int)(ceny-size);
    up_y = (int)(ceny+size);

    for (i=low_y; i<up_y; i++) {
	for (j=low_x; j<up_x; j++) {
	    rsq = (cenx - j)*(cenx - j) + (ceny - i)*(ceny - i);
	    gauss = a*expf(-1.0*rsq/(sigma*sigma));
	    if (gauss > 255)
		gauss = 255;
	    buffer[i*naxes[0]+j] += (char)gauss;
	}
    }

    return 1;
}

int main() {

    dc1394camera_t *camera;
    char *buffer, *buffer2, *average;
    fitsfile *fptr;
    int i, j, f, status, nimages, anynul, nboxes, test;
    char fitsfile[256], xpastr[256];
    char *froot, *timestr;
    FILE *init, *out;
    float xx = 0.0, yy = 0.0, xsum = 0.0, ysum = 0.0;
    double dist[2000], sig[2000], dist_l[2000], sig_l[2000];
    double mean, var, var_l, avesig;
    double seeing_short, seeing_long, seeing_ave;
    struct timeval start_time, end_time;
    time_t start_sec, end_sec;
    suseconds_t start_usec, end_usec;
    float elapsed_time, fps;

    dc1394_t * dc;
    dc1394camera_list_t * list;
    dc1394error_t err;

    unsigned int min_bytes, max_bytes, max_height, max_width;
    unsigned int actual_bytes, winleft, wintop;
    uint64_t total_bytes = 0;
    char *names[NXPA];
    char *messages[NXPA];

    double d = 0.060;
    double r = 0.130;

    XPA xpa;
    xpa = XPAOpen(NULL);

    stderr = freopen("measure_seeing.log", "w", stderr);

    status = 0;
    anynul = 0;
    naxes[0] = 320;
    naxes[1] = 240;
    fpixel = 1;

    nelements = naxes[0]*naxes[1];

    dc = dc1394_new();
    if (!dc)
	return 1;
    err = dc1394_camera_enumerate(dc, &list);
    DC1394_ERR_RTN(err, "Failed to enumerate cameras.");

    if (list->num == 0) {
	dc1394_log_error("No cameras found.");
	return 1;
    }

    camera = dc1394_camera_new(dc, list->ids[0].guid);
    if (!camera) {
	dc1394_log_error("Failed to initialize camera with guid %"PRIx64".",
			 list->ids[0].guid);
	return 1;
    }
    dc1394_camera_free_list(list);

    printf("Using camera with GUID %"PRIx64"\n", camera->guid);

    // need to use legacy firewire400 mode for now.  800 not quite reliable.
    dc1394_video_set_iso_speed(camera, DC1394_ISO_SPEED_400);

    // configure camera for format7
    err = dc1394_video_set_mode(camera, DC1394_VIDEO_MODE_FORMAT7_1);
    DC1394_ERR_CLN_RTN(err, dc1394_camera_free(camera), "Can't choose format7_0");
    printf("I: video mode is format7_0\n");

    err = dc1394_format7_get_max_image_size(camera, DC1394_VIDEO_MODE_FORMAT7_1,
					    &max_width, &max_height);
    DC1394_ERR_CLN_RTN(err,dc1394_camera_free (camera),"cannot get max image size for format7_0");
    printf ("I: max image size is: height = %d, width = %d\n", max_height, max_width);
    printf ("I: current image size is: height = %ld, width = %ld\n", naxes[1], naxes[0]);

    //winleft = (max_width - naxes[0])/2;
    //wintop = (max_height - naxes[1])/2;
    winleft = 0;
    wintop = 0;

    err = dc1394_format7_set_roi(camera,
				 DC1394_VIDEO_MODE_FORMAT7_1,
				 DC1394_COLOR_CODING_MONO8,
				 DC1394_USE_MAX_AVAIL,
				 winleft, wintop, // left, top
				 naxes[0], naxes[1]);
    DC1394_ERR_CLN_RTN(err, dc1394_camera_free(camera), "Can't set ROI.");
    printf("I: ROI is (%d, %d) - (%ld, %ld)\n",
	   winleft, wintop, winleft+naxes[0], wintop+naxes[1]);

    err = dc1394_format7_get_total_bytes(camera, DC1394_VIDEO_MODE_FORMAT7_1, &total_bytes);
    DC1394_ERR_CLN_RTN(err, dc1394_camera_free(camera), "Can't get total bytes.");
    printf("I: total bytes per frame are %"PRIu64"\n", total_bytes);

    err = dc1394_capture_setup(camera, 16, DC1394_CAPTURE_FLAGS_DEFAULT);
    DC1394_ERR_CLN_RTN(err, dc1394_camera_free(camera), "Error capturing.");

    // start the camera up
    err = dc1394_video_set_transmission(camera, DC1394_ON);
    if (err != DC1394_SUCCESS) {
	dc1394_log_error("Unable to start camera iso transmission.");
	dc1394_capture_stop(camera);
	dc1394_camera_free(camera);
	exit(1);
    }

    out = fopen("seeing.dat", "a");
    init = fopen("init_cen_all", "r");
    i = 0;
    while (fscanf(init, "%f %f\n", &xx, &yy) != EOF) {
	box[i].x = xx;
	box[i].cenx = xx;
	box[i].y = yy;
	box[i].ceny = yy;
	box[i].r = boxsize/2.0;
	i++;
    }
    nboxes = i;
    fclose(init);

    back.r = 80;
    back.width = 10;

    /* allocate the buffers */
    if (!(buffer = malloc(nelements*sizeof(char)))) {
	printf("Couldn't Allocate Image Buffer\n");
	exit(-1);
    }
    if (!(buffer2 = malloc(nelements*sizeof(char)))) {
	printf("Couldn't Allocate 2nd Image Buffer\n");
	exit(-1);
    }
    if (!(average = malloc(nelements*sizeof(char)))) {
	printf("Couldn't Allocate Average Image Buffer\n");
	exit(-1);
    }

    nimages = 2000;
    froot = "seeing.fits";

    gettimeofday(&start_time, NULL);

    for (f=0; f<nimages; f++) {
	/* first do a single exposure */
	grab_frame(camera, buffer, nelements*sizeof(char));
	// add_gaussian(buffer, 175.0, 130.0, 50.0, 3.0);
	// add_gaussian(buffer, 155.0, 115.0, 50.0, 3.0);

	// find center of star images and calculate background
	xsum = 0.0;
	ysum = 0.0;
	for (i=0; i<nboxes; i++) {
	    xsum += box[i].cenx;
	    ysum += box[i].ceny;
	}
	back.x = xsum/nboxes;
	back.y = ysum/nboxes;
	background(buffer, naxes[0], naxes[1]);

	for (i=0; i<nboxes; i++) {
	    box[i].noise = back.sigma;
	    box[i].r = boxsize/2.0;
	    //centroid(buffer, i);
	    //box[i].r = boxsize/3.0;
	    //centroid(buffer, i);
	    //box[i].r = boxsize/4.0;
	    centroid(buffer, naxes[0], naxes[1], i);
	}

	dist[f] = stardist(0, 1);
	sig[f] = box[0].sigmaxy*box[0].sigmaxy + box[1].sigmaxy*box[1].sigmaxy;

	/* now average two exposures */
	grab_frame(camera, buffer, nelements*sizeof(char));
	grab_frame(camera, buffer2, nelements*sizeof(char));

	for (j=0; j<nelements; j++) {
	    test = buffer[j]+buffer2[j];
	    if (test <= 254) {
		average[j] = buffer[j]+buffer2[j];
	    } else {
		average[j] = 254;
	    }
	}
	// add_gaussian(average, 175.0, 130.0, 50.0, 3.0);
	// add_gaussian(average, 155.0, 115.0, 50.0, 3.0);

	xsum = 0.0;
	ysum = 0.0;
	for (i=0; i<nboxes; i++) {
	    xsum += box[i].cenx;
	    ysum += box[i].ceny;
	}
	back.x = xsum/nboxes;
	back.y = ysum/nboxes;
	background(average, naxes[0], naxes[1]);

	for (i=0; i<nboxes; i++) {
	    box[i].noise = back.sigma;
	    box[i].r = boxsize/2.0;
	    //centroid(average, i);
	    //box[i].r = boxsize/3.0;
	    //centroid(average, i);
	    //box[i].r = boxsize/4.0;
	    centroid(average, naxes[0], naxes[1], i);
	}

	dist_l[f] = stardist(0, 1);
	sig_l[f] = box[0].sigmaxy*box[0].sigmaxy + box[1].sigmaxy*box[1].sigmaxy;

	if (f % 40 == 0) {
	    status = XPASet(xpa, "ds9", "array [xdim=320,ydim=240,bitpix=8]", "ack=false",
			    average, nelements, names, messages, NXPA);
	    sprintf(xpastr, "image; box %f %f %d %d 0.0",
		    box[0].x, box[0].y, boxsize, boxsize);
	    status = XPASet(xpa, "ds9", "regions", "ack=false",
			    xpastr, strlen(xpastr), names, messages, NXPA);
	    sprintf(xpastr, "image; box %f %f %d %d 0.0",
		    box[1].x, box[1].y, boxsize, boxsize);
	    status = XPASet(xpa, "ds9", "regions", "ack=false",
			    xpastr, strlen(xpastr), names, messages, NXPA);
	}

    }

    gettimeofday(&end_time, NULL);
    printf("End capture.\n");

    /*-----------------------------------------------------------------------
     *  stop data transmission
     *-----------------------------------------------------------------------*/
    start_sec = start_time.tv_sec;
    start_usec = start_time.tv_usec;
    end_sec = end_time.tv_sec;
    end_usec = end_time.tv_usec;

    elapsed_time = (float)((end_sec + 1.0e-6*end_usec) - (start_sec + 1.0e-6*start_usec));
    fps = 3*nimages/elapsed_time;
    printf("Elapsed time = %g seconds.\n", elapsed_time);
    printf("Framerate = %g fps.\n", fps);

    err=dc1394_video_set_transmission(camera,DC1394_OFF);
    DC1394_ERR_RTN(err,"couldn't stop the camera?");

    /*
    sprintf(fitsfile, "!%s", froot);
    fits_create_file(&fptr, fitsfile, &status);
    fits_create_img(fptr, BYTE_IMG, 2, naxes, &status);
    fits_write_img(fptr, TBYTE, fpixel, nelements, buffer, &status);
    fits_close_file(fptr, &status);
    fits_report_error(stderr, status);
    */

    /* analyze short exposure */
    printf("\t SHORT EXPOSURE\n");
    mean = gsl_stats_mean(dist, 1, nimages);
    avesig = gsl_stats_mean(sig, 1, nimages);
    printf("mean = %f, avesig = %f\n", mean, avesig);

    printf("\n");

    var = gsl_stats_variance_m(dist, 1, nimages, mean);
    var = var - avesig;
    seeing_short = seeing(var, d, r);
    printf("sigma = %f, seeing = %f\n", sqrt(var), seeing_short);

    /* analyze long exposure */
    printf("\t LONG EXPOSURE\n");
    mean = gsl_stats_mean(dist_l, 1, nimages);
    avesig = gsl_stats_mean(sig_l, 1, nimages);
    printf("mean_l = %f, avesig_l = %f\n", mean, avesig);

    printf("\n");

    var_l = gsl_stats_variance_m(dist_l, 1, nimages, mean);
    var_l = var_l - avesig;
    seeing_long = seeing(var_l, d, r);
    printf("sigma_l = %f, seeing_l = %f\n", sqrt(var_l), seeing_long);

    seeing_ave = pow(seeing_short, 1.75)*pow(seeing_long,-0.75);
    printf("Exposure corrected seeing = %4.2f\"\n\n", seeing_ave);

    timestr = ctime(&end_sec);
    fprintf(out, "%s %f %f %f %f %f\n", timestr, var, var_l, seeing_short, seeing_long, seeing_ave);

    init = fopen("init_cen_all", "w");
    for (i=0; i<nboxes; i++) {
	fprintf(init, "%f %f\n", box[i].cenx, box[i].ceny);
    }
    fclose(init);

    fclose(out);

    return (status);
}

max_stars_repo_name: marissakotze/timDIMM

max_issues_repo_name: marissakotze/timDIMM

max_forks_repo_name: marissakotze/timDIMM


text:
#ifndef QUAC_P_H_
#define QUAC_P_H_
#include <petsc.h>
extern int  petsc_initialized;
PetscLogEvent add_lin_event,add_to_ham_event,add_lin_recovery_event,add_encoded_gate_to_circuit_event;
PetscLogEvent _qc_event_function_event,_qc_postevent_function_event,_apply_gate_event;
PetscClassId quac_class_id;
PetscLogStage pre_solve_stage,solve_stage,post_solve_stage;
#endif

max_stars_repo_name: sgulania/QuaC

max_issues_repo_name: sgulania/QuaC

max_forks_repo_name: sgulania/QuaC


text:
//**************************************************************************\
//* This file is property of and copyright by the ALICE Project            *\
//* ALICE Experiment at CERN, All rights reserved.                         *\
//*                                                                        *\
//* Primary Authors: Matthias Richter <Matthias.Richter@ift.uib.no>        *\
//*                  for The ALICE HLT Project.                            *\
//*                                                                        *\
//* Permission to use, copy, modify and distribute this software and its   *\
//* documentation strictly for non-commercial purposes is hereby granted   *\
//* without fee, provided that the above copyright notice appears in all   *\
//* copies and that both the copyright notice and this permission notice   *\
//* appear in the supporting documentation. The authors make no claims     *\
//* about the suitability of this software for any purpose. It is          *\
//* provided "as is" without express or implied warranty.                  *\
//**************************************************************************

/// \file GPUO2InterfaceRefit.h
/// \author David Rohr

#ifndef GPUO2INTERFACEREFIT_H
#define GPUO2INTERFACEREFIT_H

// Some defines denoting that we are compiling for O2
#ifndef HAVE_O2HEADERS
#define HAVE_O2HEADERS
#endif
#ifndef GPUCA_TPC_GEOMETRY_O2
#define GPUCA_TPC_GEOMETRY_O2
#endif
#ifndef GPUCA_O2_INTERFACE
#define GPUCA_O2_INTERFACE
#endif

#include <memory>
#include <vector>
#include <gsl/span>

namespace o2::base
{
template <typename value_T>
class PropagatorImpl;
using Propagator = PropagatorImpl<float>;
} // namespace o2::base
namespace o2::dataformats
{
template <typename FirstEntry, typename NElem>
class RangeReference;
}
namespace o2::tpc
{
using TPCClRefElem = uint32_t;
using TrackTPCClusRef = o2::dataformats::RangeReference<uint32_t, uint16_t>;
class TrackTPC;
struct ClusterNativeAccess;
} // namespace o2::tpc
namespace o2::track
{
template <typename value_T>
class TrackParametrizationWithError;
using TrackParCovF = TrackParametrizationWithError<float>;
using TrackParCov = TrackParCovF;
} // namespace o2::track

namespace o2::gpu
{
class GPUParam;
class GPUTrackingRefit;
class TPCFastTransform;
class GPUO2InterfaceRefit
{
 public:
  // Must initialize with:
  // - In any case: Cluster Native access structure (cl), TPC Fast Transformation instance (trans), solenoid field (bz), TPC Track hit references (trackRef)
  // - Either the shared cluster map (sharedmap) or the vector of tpc tracks (trks) to build the shared cluster map internally
  // - o2::base::Propagator (p) in case RefitTrackAsTrackParCov is to be used

  GPUO2InterfaceRefit(const o2::tpc::ClusterNativeAccess* cl, const TPCFastTransform* trans, float bz, const o2::tpc::TPCClRefElem* trackRef, const unsigned char* sharedmap = nullptr, const std::vector<o2::tpc::TrackTPC>* trks = nullptr, o2::base::Propagator* p = nullptr);
  ~GPUO2InterfaceRefit();

  int RefitTrackAsGPU(o2::tpc::TrackTPC& trk, bool outward = false, bool resetCov = false);
  int RefitTrackAsTrackParCov(o2::tpc::TrackTPC& trk, bool outward = false, bool resetCov = false);
  int RefitTrackAsGPU(o2::track::TrackParCov& trk, const o2::tpc::TrackTPCClusRef& clusRef, float time0, float* chi2 = nullptr, bool outward = false, bool resetCov = false);
  int RefitTrackAsTrackParCov(o2::track::TrackParCov& trk, const o2::tpc::TrackTPCClusRef& clusRef, float time0, float* chi2 = nullptr, bool outward = false, bool resetCov = false);
  void setGPUTrackFitInProjections(bool v = true);
  void setTrackReferenceX(float v);
  void setIgnoreErrorsAtTrackEnds(bool v);

  static void fillSharedClustersMap(const o2::tpc::ClusterNativeAccess* cl, const gsl::span<const o2::tpc::TrackTPC> trks, const o2::tpc::TPCClRefElem* trackRef, unsigned char* shmap);

 private:
  std::unique_ptr<GPUTrackingRefit> mRefit;
  std::unique_ptr<GPUParam> mParam;
  std::vector<unsigned char> mSharedMap;
};
} // namespace o2::gpu

#endif


text:
//
//     smctc.h
//
//     The principle header file for the SMC template class.
//
//     Copyright Adam Johansen, 2008
//
//
//   This file is part of SMCTC.
//
//   SMCTC is free software: you can redistribute it and/or modify
//   it under the terms of the GNU General Public License as published by
//   the Free Software Foundation, either version 3 of the License, or
//   (at your option) any later version.
//
//   SMCTC is distributed in the hope that it will be useful,
//   but WITHOUT ANY WARRANTY; without even the implied warranty of
//   MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
//   GNU General Public License for more details.
//
//   You should have received a copy of the GNU General Public License
//   along with SMCTC.  If not, see <http://www.gnu.org/licenses/>.
//

/// \mainpage smctc -- A Template Class Library for SMC Simulation
///
/// \version 1.0
/// \author Adam M. Johansen
///
/// \section intro_sec Summary
///
/// The SMC template class library (SMCTC) is intended to be used to implement SMC
/// algorithms ranging from simple particle filter to complicated SMC algorithms
/// from simple particle filters to the SMC samplers of Del Moral, Doucet and Jasra (2006)
/// within a generic framework.
///
/// \section license License
///
/// SMCTC is free software: you can redistribute it and/or modify
/// it under the terms of the GNU General Public License as published by
/// the Free Software Foundation, either version 3 of the License, or
/// (at your option) any later version.
///
/// SMCTC is distributed in the hope that it will be useful,
/// but WITHOUT ANY WARRANTY; without even the implied warranty of
/// MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
/// GNU General Public License for more details.
///
/// You should have received a copy of the GNU General Public License
/// along with SMCTC.  If not, see <http://www.gnu.org/licenses/>.
///
/// The software is still in development but is thought to be sufficiently fully-featured and
/// stable for use.
///
/// \section using_sec Using the Template Library
///
/// In order to use the template libary it is necessary to place the set of header files in
/// a suitable include directory and then to instantiate the various classes with the
/// appropriate type for your sampler.
///
/// For implementation details, please see the accompanying user guide.
///
///
/// \section example_sec Examples
///
/// \subsection A Simple Particle Filter
///
/// This example provides a simple particle filter for the approximate filtering of
/// nonlinear, nongaussian state space models.
///
/// \subsection SMC Samplers for Rare Event Simulation
///
/// The main.cc file provides an example of a use of the template library to estimate
/// Gaussian tail probabilities within an SMC samplers framework by making use of a sequence
/// of intermediate distributions defined by introducing a potential function which
/// varies from flat to being very close to an indicator function on the tail event.
///
/// This example is taken from from Johansen, Del Moral and Doucet (2006).
///
/// \subsection Acknowledgements
///
/// Thanks are due to Edmund Jackson and Mark Briers for testing SMCTC on a variety of platforms.
/// The Visual C++ project and solution files are kindly provided by Mark Briers.
///

//! \file
//! \brief The main header file for SMCTC.
//!
//! This file serves as an interface between user-space programs and the SMCTC library.
//! This is the only header file which applications need to include to make use of the library (it includes
//! such additional files as are necessary).

#ifndef __SMC_TDSMC_HH

#define __SMC_TDSMC_HH 1.0

#include <cmath>
#include <cstdlib>
#include <iostream>

#include <gsl/gsl_rng.h>

#include "smc-exception.h"
#include "sampler.h"

/// The Sequential Monte Carlo namespace

///
///    The classes and functions within this namespace are intended to be used for producing
///    implemenetations of SMC samplers and related simulation techniques.

namespace smc {}

/// The standard namespace

///     The classes provided within the standard libraries reside within this namespace and
///     the TDSMC class library adds a number of additional operator overloads to some of
///     the standard classes to allow them to deal with our structures.

namespace std {}
#endif


text:
/**
 *
 * @file qwrapper_slantr.c
 *
 *  PLASMA core_blas quark wrapper
 *  PLASMA is a software package provided by Univ. of Tennessee,
 *  Univ. of California Berkeley and Univ. of Colorado Denver
 *
 * @version 2.6.0
 * @author Mathieu Faverge
 * @date 2010-11-15
 * @generated s Tue Jan  7 11:44:57 2014
 *
 **/
#include <lapacke.h>
#include "common.h"

void
CORE_slantr_quark(Quark *quark);
void
CORE_slantr_f1_quark(Quark *quark);

/***************************************************************************//**
 *
 **/
void QUARK_CORE_slantr(Quark *quark, Quark_Task_Flags *task_flags,
                       PLASMA_enum norm, PLASMA_enum uplo, PLASMA_enum diag, int M, int N,
                       const float *A, int LDA, int szeA,
                       int szeW, float *result)
{
    szeW = max(1, szeW);
    DAG_CORE_LANGE;
    QUARK_Insert_Task(quark, CORE_slantr_quark, task_flags,
        sizeof(PLASMA_enum),                &norm,  VALUE,
        sizeof(PLASMA_enum),                &uplo,  VALUE,
        sizeof(PLASMA_enum),                &diag,  VALUE,
        sizeof(int),                        &M,     VALUE,
        sizeof(int),                        &N,     VALUE,
        sizeof(float)*szeA,     A,             INPUT,
        sizeof(int),                        &LDA,   VALUE,
        sizeof(float)*szeW,                 NULL,          SCRATCH,
        sizeof(float),                      result,        OUTPUT,
        0);
}

/***************************************************************************//**
 *
 **/
void QUARK_CORE_slantr_f1(Quark *quark, Quark_Task_Flags *task_flags,
                          PLASMA_enum norm, PLASMA_enum uplo, PLASMA_enum diag, int M, int N,
                          const float *A, int LDA, int szeA,
                          int szeW, float *result,
                          float *fake, int szeF)
{
    szeW = max(1, szeW);
    DAG_CORE_LANGE;

    if ( result == fake ) {
        QUARK_Insert_Task(quark, CORE_slantr_quark, task_flags,
            sizeof(PLASMA_enum),                &norm,  VALUE,
            sizeof(PLASMA_enum),                &uplo,  VALUE,
            sizeof(PLASMA_enum),                &diag,  VALUE,
            sizeof(int),                        &M,     VALUE,
            sizeof(int),                        &N,     VALUE,
            sizeof(float)*szeA,     A,             INPUT,
            sizeof(int),                        &LDA,   VALUE,
            sizeof(float)*szeW,                 NULL,          SCRATCH,
            sizeof(float),                      result,        OUTPUT | GATHERV,
            0);
    } else {
        QUARK_Insert_Task(quark, CORE_slantr_f1_quark, task_flags,
            sizeof(PLASMA_enum),                &norm,  VALUE,
            sizeof(PLASMA_enum),                &uplo,  VALUE,
            sizeof(PLASMA_enum),                &diag,  VALUE,
            sizeof(int),                        &M,     VALUE,
            sizeof(int),                        &N,     VALUE,
            sizeof(float)*szeA,     A,             INPUT,
            sizeof(int),                        &LDA,   VALUE,
            sizeof(float)*szeW,                 NULL,          SCRATCH,
            sizeof(float),                      result,        OUTPUT,
            sizeof(float)*szeF,                 fake,          OUTPUT | GATHERV,
            0);
    }
}

/***************************************************************************//**
 *
 **/
#if defined(PLASMA_HAVE_WEAK)
#pragma weak CORE_slantr_quark = PCORE_slantr_quark
#define CORE_slantr_quark PCORE_slantr_quark
#endif
void CORE_slantr_quark(Quark *quark)
{
    float *normA;
    PLASMA_enum norm, uplo, diag;
    int M;
    int N;
    float *A;
    int LDA;
    float *work;

    quark_unpack_args_9(quark, norm, uplo, diag, M, N, A, LDA, work, normA);
    *normA = LAPACKE_slantr_work(
        LAPACK_COL_MAJOR,
        lapack_const(norm),
        lapack_const(uplo),
        lapack_const(diag),
        M, N, A, LDA, work);
}

/***************************************************************************//**
 *
 **/
#if defined(PLASMA_HAVE_WEAK)
#pragma weak CORE_slantr_f1_quark = PCORE_slantr_f1_quark
#define CORE_slantr_f1_quark PCORE_slantr_f1_quark
#endif
void CORE_slantr_f1_quark(Quark *quark)
{
    float *normA;
    PLASMA_enum norm, uplo, diag;
    int M;
    int N;
    float *A;
    int LDA;
    float *work;
    float *fake;

    quark_unpack_args_10(quark, norm, uplo, diag, M, N, A, LDA, work, normA, fake);
    *normA = LAPACKE_slantr_work(
        LAPACK_COL_MAJOR,
        lapack_const(norm),
        lapack_const(uplo),
        lapack_const(diag),
        M, N, A, LDA, work);
}


text:
/* gsl_spmatrix.h
 *
 * Copyright (C) 2012-2014 Patrick Alken
 *
 * This program is free software; you can redistribute it and/or modify
 * it under the terms of the GNU General Public License as published by
 * the Free Software Foundation; either version 3 of the License, or (at
 * your option) any later version.
 *
 * This program is distributed in the hope that it will be useful, but
 * WITHOUT ANY WARRANTY; without even the implied warranty of
 * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the GNU
 * General Public License for more details.
 *
 * You should have received a copy of the GNU General Public License
 * along with this program; if not, write to the Free Software
 * Foundation, Inc., 51 Franklin Street, Fifth Floor, Boston, MA 02110-1301, USA.
 */

#ifndef __GSL_SPMATRIX_H__
#define __GSL_SPMATRIX_H__

#include <stdlib.h>

#include <gsl/gsl_math.h>
#include <gsl/gsl_vector.h>
#include <gsl/gsl_matrix.h>

#undef __BEGIN_DECLS
#undef __END_DECLS
#ifdef __cplusplus
# define __BEGIN_DECLS extern "C" {
# define __END_DECLS }
#else
# define __BEGIN_DECLS /* empty */
# define __END_DECLS /* empty */
#endif

__BEGIN_DECLS

/*
 * Binary tree data structure for storing sparse matrix elements
 * in triplet format. This is used for efficiently detecting
 * duplicates and element retrieval via gsl_spmatrix_get
 */
typedef struct
{
  void *tree;       /* tree structure */
  void *node_array; /* preallocated array of tree nodes */
  size_t n;         /* number of tree nodes in use (<= nzmax) */
} gsl_spmatrix_tree;

/*
 * Triplet format:
 *
 * If data[n] = A_{ij}, then:
 *   i = A->i[n]
 *   j = A->p[n]
 *
 * Compressed column format (CCS):
 *
 * If data[n] = A_{ij}, then:
 *   i = A->i[n]
 *   A->p[j] <= n < A->p[j+1]
 * so that column j is stored in
 * [ data[p[j]], data[p[j] + 1], ..., data[p[j+1] - 1] ]
 *
 * Compressed row format (CRS):
 *
 * If data[n] = A_{ij}, then:
 *   j = A->i[n]
 *   A->p[i] <= n < A->p[i+1]
 * so that row i is stored in
 * [ data[p[i]], data[p[i] + 1], ..., data[p[i+1] - 1] ]
 */

typedef struct
{
  size_t size1;  /* number of rows */
  size_t size2;  /* number of columns */

  /* i (size nzmax) contains:
   *
   * Triplet/CCS: row indices
   * CRS: column indices
   */
  size_t *i;

  double *data;  /* matrix elements of size nzmax */

  /*
   * p contains the column indices (triplet) or column pointers (compcol)
   *
   * triplet: p[n] = column number of element data[n]
   * CCS:     p[j] = index in data of first non-zero element in column j
   * CRS:     p[i] = index in data of first non-zero element in row i
   */
  size_t *p;

  size_t nzmax;  /* maximum number of matrix elements */
  size_t nz;     /* number of non-zero values in matrix */

  gsl_spmatrix_tree *tree_data; /* binary tree for sorting triplet data */

  /*
   * workspace of size MAX(size1,size2)*MAX(sizeof(double),sizeof(size_t))
   * used in various routines
   */
  void *work;

  size_t sptype; /* sparse storage type */
} gsl_spmatrix;

#define GSL_SPMATRIX_TRIPLET      (0)
#define GSL_SPMATRIX_CCS          (1)
#define GSL_SPMATRIX_CRS          (2)

#define GSL_SPMATRIX_ISTRIPLET(m) ((m)->sptype == GSL_SPMATRIX_TRIPLET)
#define GSL_SPMATRIX_ISCCS(m)     ((m)->sptype == GSL_SPMATRIX_CCS)
#define GSL_SPMATRIX_ISCRS(m)     ((m)->sptype == GSL_SPMATRIX_CRS)

/*
 * Prototypes
 */

gsl_spmatrix *gsl_spmatrix_alloc(const size_t n1, const size_t n2);
gsl_spmatrix *gsl_spmatrix_alloc_nzmax(const size_t n1, const size_t n2,
                                       const size_t nzmax, const size_t flags);
void gsl_spmatrix_free(gsl_spmatrix *m);
int gsl_spmatrix_realloc(const size_t nzmax, gsl_spmatrix *m);
int gsl_spmatrix_set_zero(gsl_spmatrix *m);
size_t gsl_spmatrix_nnz(const gsl_spmatrix *m);
int gsl_spmatrix_compare_idx(const size_t ia, const size_t ja,
                             const size_t ib, const size_t jb);
int gsl_spmatrix_tree_rebuild(gsl_spmatrix * m);

/* spcopy.c */
int gsl_spmatrix_memcpy(gsl_spmatrix *dest, const gsl_spmatrix *src);

/* spgetset.c */
double gsl_spmatrix_get(const gsl_spmatrix *m, const size_t i,
                        const size_t j);
int gsl_spmatrix_set(gsl_spmatrix *m, const size_t i, const size_t j,
                     const double x);
double *gsl_spmatrix_ptr(gsl_spmatrix *m, const size_t i, const size_t j);

/* spcompress.c */
gsl_spmatrix *gsl_spmatrix_compcol(const gsl_spmatrix *T);
gsl_spmatrix *gsl_spmatrix_ccs(const gsl_spmatrix *T);
gsl_spmatrix *gsl_spmatrix_crs(const gsl_spmatrix *T);
void gsl_spmatrix_cumsum(const size_t n, size_t *c);

/* spio.c */
int gsl_spmatrix_fprintf(FILE *stream, const gsl_spmatrix *m,
                         const char *format);
gsl_spmatrix * gsl_spmatrix_fscanf(FILE *stream);
int gsl_spmatrix_fwrite(FILE *stream, const gsl_spmatrix *m);
int gsl_spmatrix_fread(FILE *stream, gsl_spmatrix *m);

/* spoper.c */
int gsl_spmatrix_scale(gsl_spmatrix *m, const double x);
int gsl_spmatrix_minmax(const gsl_spmatrix *m, double *min_out,
                        double *max_out);
int gsl_spmatrix_add(gsl_spmatrix *c, const gsl_spmatrix *a,
                     const gsl_spmatrix *b);
int gsl_spmatrix_d2sp(gsl_spmatrix *S, const gsl_matrix *A);
int gsl_spmatrix_sp2d(gsl_matrix *A, const gsl_spmatrix *S);

/* spprop.c */
int gsl_spmatrix_equal(const gsl_spmatrix *a, const gsl_spmatrix *b);

/* spswap.c */
int gsl_spmatrix_transpose(gsl_spmatrix * m);
int gsl_spmatrix_transpose2(gsl_spmatrix * m);
int gsl_spmatrix_transpose_memcpy(gsl_spmatrix *dest, const gsl_spmatrix *src);

__END_DECLS

#endif /* __GSL_SPMATRIX_H__ */

max_stars_repo_name: Tjoppen/fmigo

max_issues_repo_name: Brian-ning/HMNE

max_forks_repo_name: Brian-ning/HMNE


text:
/**
 *
 * @file example_dpotrf.c
 *
 *  PLASMA testing routines
 *  PLASMA is a software package provided by Univ. of Tennessee,
 *  Univ. of California Berkeley and Univ. of Colorado Denver
 *
 * @brief Example of Cholesky factorization
 *
 * @version 2.6.0
 * @author Bilel Hadri
 * @date 2010-11-15
 * @generated d Tue Jan  7 11:45:20 2014
 *
 **/
#include <stdlib.h>
#include <stdio.h>
#include <string.h>
#include <math.h>
#include <time.h>

#include <plasma.h>
#include <cblas.h>
#include <lapacke.h>
#include <core_blas.h>
#include "testing_dmain.h"

static int check_factorization(int, double*, double*, int, int);

//int IONE=1;
//int ISEED[4] = {0,0,0,1};   /* initial seed for dlarnv() */

static void GENMAT_SYM_FULL(int m, double *A)
{
	srand48(time(NULL));

	int j;
	for (j = 0; j < m; ++j ) {
		int i;
		for( i = j; i < m; ++i ) {
			double dran = drand48();
			A[j*m+i] = A[i*m+j] = dran;
		}
  	}
	for(j = 0; j < m; ++j)
		A[j*m+j] += 10 * m;
}

int testing_dpotrf(int argc, char **argv)
{

	if (argc != 2){
		fprintf(stderr,"POTRF: N LDA\n");
		return -1;
	}

    int N     = atoi(argv[0]);
    int LDA   = atoi(argv[1]);
    int info_factorization;

	PLASMA_Set(PLASMA_TILE_SIZE, 388);
    double *A1   = (double *)malloc(LDA*N*sizeof(double));
	#pragma omp register ([LDA*N]A1)
    double *A2   = (double *)malloc(LDA*N*sizeof(double));
	#pragma omp register ([LDA*N]A2)

    /* Check if unable to allocate memory */
    if ((!A1)||(!A2)){
        printf("Out of Memory \n ");
        return 0;
    }

    /* Initialize A1 and A2 for Symmetric Positive Matrix */
	GENMAT_SYM_FULL(N, A1);
	int i;
	for(i = 0; i < N*LDA; ++i){
		A2[i] = A1[i];
	}

    /* Plasma routines */
    PLASMA_dpotrf(PlasmaUpper, N, A2, LDA);

    /* Check the factorization */
    info_factorization = check_factorization( N, A1, A2, LDA, PlasmaUpper);

    if ( info_factorization != 0 )
       printf("-- Error in DPOTRF example ! \n");
    else
       printf("-- Run of DPOTRF example successful ! \n");

    free(A1); free(A2);

    return 0;
}

static int check_factorization(int N, double *A1, double *A2, int LDA, int uplo)
{
    double Anorm, Rnorm;
    double alpha;
    int info_factorization;
    int i,j;
    double eps;

    eps = LAPACKE_dlamch_work('e');

    double *Residual = (double *)malloc(N*N*sizeof(double));
    double *L1       = (double *)malloc(N*N*sizeof(double));
    double *L2       = (double *)malloc(N*N*sizeof(double));
    double *work              = (double *)malloc(N*sizeof(double));

    memset((void*)L1, 0, N*N*sizeof(double));
    memset((void*)L2, 0, N*N*sizeof(double));

    alpha= 1.0;

    LAPACKE_dlacpy_work(LAPACK_COL_MAJOR,' ', N, N, A1, LDA, Residual, N);

    /* Dealing with L'L or U'U  */
    if (uplo == PlasmaUpper){
        LAPACKE_dlacpy_work(LAPACK_COL_MAJOR,'u', N, N, A2, LDA, L1, N);
        LAPACKE_dlacpy_work(LAPACK_COL_MAJOR,'u', N, N, A2, LDA, L2, N);
        cblas_dtrmm(CblasColMajor, CblasLeft, CblasUpper, CblasTrans, CblasNonUnit, N, N, (alpha), L1, N, L2, N);
    }
    else{
        LAPACKE_dlacpy_work(LAPACK_COL_MAJOR,'l', N, N, A2, LDA, L1, N);
        LAPACKE_dlacpy_work(LAPACK_COL_MAJOR,'l', N, N, A2, LDA, L2, N);
        cblas_dtrmm(CblasColMajor, CblasRight, CblasLower, CblasTrans, CblasNonUnit, N, N, (alpha), L1, N, L2, N);
    }

    /* Compute the Residual || A -L'L|| */
    for (i = 0; i < N; i++)
        for (j = 0; j < N; j++)
           Residual[j*N+i] = L2[j*N+i] - Residual[j*N+i];

    Rnorm = LAPACKE_dlange_work(LAPACK_COL_MAJOR, lapack_const(PlasmaInfNorm), N, N, Residual, N, work);
    Anorm = LAPACKE_dlange_work(LAPACK_COL_MAJOR, lapack_const(PlasmaInfNorm), N, N, A1, LDA, work);

    printf("============\n");
    printf("Checking the Cholesky Factorization \n");
    printf("-- ||L'L-A||_oo/(||A||_oo.N.eps) = %e \n",Rnorm/(Anorm*N*eps));

    if ( isnan(Rnorm/(Anorm*N*eps)) || (Rnorm/(Anorm*N*eps) > 10.0) ){
        printf("-- Factorization is suspicious ! \n");
        info_factorization = 1;
    }
    else{
        printf("-- Factorization is CORRECT ! \n");
        info_factorization = 0;
    }

    free(Residual); free(L1); free(L2); free(work);

    return info_factorization;
}

int timing_dpotrf(int argc, char **argv)
{

    int N;
    int LDA;
	int bs;
	int rep;

	double start, end;
	double elapsed;
	FILE *log;
	int num_threads;

	if (argc != 4 ) {
		fprintf(stderr,"DPOTRF n lda bs rep\n");
		exit(1);
	}

	sscanf(argv[0], "%d", &N);
	sscanf(argv[1], "%d", &LDA);
	sscanf(argv[2], "%d", &bs);
	sscanf(argv[3], "%d", &rep);

    double *A = malloc(LDA*N*sizeof(double));
	#pragma omp register ([LDA*N]A)

    /* Check if unable to allocate memory */
    if (!A){
        printf("Out of Memory \n ");
        return 0;
    }

    /* Plasma Initialize */
	PLASMA_Set( PLASMA_TILE_SIZE, bs);

	GENMAT_SYM_FULL(N, A);

	elapsed = 0.0;
	int i;
	for ( i = 0; i < rep; i++ ) {
		start = gtime();
		PLASMA_dpotrf(PlasmaUpper, N, A, LDA);
		end = gtime();
		elapsed += end - start;
	}

	num_threads = omp_get_max_threads();
	dump_info("plasma_dpotrf.log", num_threads, elapsed, rep);

    free(A);

    return 0;

}