text:
/*
 * HelloWorld.cpp
 *
 *  Created on: 22 August 2016
 *      Author: Dr. Corneliu Arsene
 *
 *   This software implements the Canonical Variates Analysis method (libhelloworld.jnilib file) for image
 *   processing: it reads the multispectral images from the same directory as where the libhelloworld.jnilib
 *   file is located  and then produces a number of improved colour images in comparison to the original
 *   multispectral images.
 *
 *
 *
 *
 *
 */

#include "jni.h"
#include <stdio.h>
#include "HelloWorld.h"

#include <math.h>
#include <assert.h>

#include <dirent.h>

#include </usr/local/Cellar/libiomp/20150701/include/libiomp/omp.h>

#include <opencv2/core/core.hpp>
#include <opencv2/highgui.hpp>

#include <opencv/cv.h>
#include <opencv/highgui.h>

#include <time.h>

#include <string.h>
#include <sys/stat.h>
#include <stdlib.h>
#include <unistd.h>
#include <string.h>

#include <iostream>

#include "gsl/matrix/gsl_matrix.h"
#include "gsl/specfunc/gsl_sf_bessel.h"
#include "gsl_math.h"
#include "gsl/gsl_blas.h"
#include <gsl/gsl_vector.h>
#include <gsl/gsl_errno.h>
#include <gsl/gsl_eigen.h>
#include <gsl/gsl_sort_vector.h>
#include <gsl/eigen/gsl_eigen.h>

using namespace cv;
using namespace std;

// * THIS FILE IS COMPILED TO A LIBRARY WHICH IS CALLED BY A JAVA LIBRARY

JNIEXPORT void JNICALL
Java_HelloWorld_print(JNIEnv *env, jobject thisObj){

    int i,j,k;
    int grp1[9];//,grouping_vector[200] ;

    char *dir;
    int depth;

    dir ="/home";
    depth=0;

    DIR * dp1,*dp3;
    struct dirent* ep,*ep3;

    char* text;
    char pattern[6];
    char* text2;
    pattern[0]='\0';
    pattern[1]='\0';
    pattern[2]='\0';
    pattern[3]='\0';
    pattern[4]='\0';
    pattern[5]='\0';

    clock_t start, end;
    double cpu_time_used;

    start = clock();
    // * DETECT HOW MANY CLASS TEXT FILES EXIST (E.G. CLASS MANUSCRIPT, CLASS UNDERWRITING, ETC)

    text="Class";

    dp1 = opendir ("./");

    int acindex,newacindex;
    acindex = 0;
    newacindex =0;

    while (ep = readdir (dp1))
    {
        text2 = ep->d_name;

        pattern[0] = text2[0];//ep->d_name[0];
        pattern[1] = text2[1];//ep->d_name[1];
        pattern[2] = text2[2];//ep->d_name[2];
        pattern[3] = text2[3];//ep->d_name[3];
        pattern[4] = text2[4];//ep->d_name[4];
        pattern[5] ='\0';

        if( strcmp(text, pattern) == 0  )
        {
            acindex = acindex +1;
        }

    }

    (void) closedir (dp1);

    // * ONCE I KNOW HOW MANY CLASSES THERE ARE THEN I MEMORIZE THE PATH TO THE RESPECTIVE CLASS FILES

    char arrayclass[acindex][2024];
    dp3 = opendir ("./");

    while (ep3 = readdir (dp3))

    {

        text2 = ep3->d_name;

        pattern[0] = text2[0];//ep->d_name[0];

        pattern[1] = text2[1];//ep->d_name[1];

        pattern[2] = text2[2];//ep->d_name[2];

        pattern[3] = text2[3];//ep->d_name[3];

        pattern[4] = text2[4];//ep->d_name[4];

        pattern[5] ='\0';

        if( strcmp(text, pattern) == 0  )
        {

            strcpy(arrayclass[newacindex],ep3->d_name);

            newacindex = newacindex +1;

        }
    }

    // * I AM PREPARING TO READ THE CLASS TEXT FILES WHICH I DETECTED ABOVE

    (void) closedir (dp3);
    char cwd[1024],cwdtemp[1024];
    getcwd(cwd, sizeof(cwd));

    strcat(cwd,"/");

    FILE *file[acindex];

    for (int j=0; j < acindex; j=j+1)
    {
        strcpy(cwdtemp,cwd);

        strcat(cwdtemp,arrayclass[j]);
        if((file[j] = fopen(cwdtemp, "r"))==NULL)       /* open a text file for reading */
            exit(1);
        printf("%s\n",cwdtemp);
    }

    int over,over1,over2,over3,over4,over5,over6;

    int xx=0;
    double dd=0;
    char str[256],str1[256],str2[256],str3[256],str4[256],str5[256],str6[256],str7[256],str8[256];

    int classoverwriting[50][7];

    double classoverwriting1[1][7];

    int classunderwriting[50][7];
    int classmanuscript[50][7];
    int classbothuo[50][7];

    int i1;

    int nrclass[acindex];

    // * I AM DETECTING HOW MANY POINTS HAS EACH CLASS, FOR EXAMPLE CLASS MANUSCRIPT MIGHT HAVE 100 POINTS WHILE CLASS OVERWRITING MIGHT HAVE 120 POINTS, ETC

    for (int j=0; j < acindex; j=j+1)
    {
        fscanf(file[j], "  %s %s %s %s %s %s \n", &str,&str1, &str2,&str3,&str4,&str5);

        i1=0;

        while (!feof(file[j]))
        {
            fscanf(file[j], "  %lf %lf %lf %lf %lf %lf %lf\n", &classoverwriting1[1][1], &classoverwriting1[1][2], &classoverwriting1[1][3], &classoverwriting1[1][4], &classoverwriting1[1][5], &classoverwriting1[1][6], &classoverwriting1[1][7]);

            i1 =i1+1;

        }
        nrclass[j]=i1;

    }

    int sumtotalnrclass;
    sumtotalnrclass = 0;

    //exit(1);

    for (int j=0; j < acindex; j=j+1)
    {
        sumtotalnrclass  = sumtotalnrclass  +  nrclass[j];
    }

    int dummyvar;

    int **classtotal = new int*[sumtotalnrclass];
    for (int i = 0; i < sumtotalnrclass; i++)
    {
        classtotal[i] = new int[7];
    }

    printf(" \n");

    double testclass1,testclass2;

    // * I AM READING THE POINTS FROM EACH CLASS FILES

    i1=0;
    for (int j=0; j < acindex; j=j+1)
    {
        fseek( file[j], 0, SEEK_SET );
        fscanf(file[j], "  %s %s %s %s %s %s \n", &str,&str1, &str2,&str3,&str4,&str5);
        //printf ("%s\n", str);
        while (!feof(file[j]))
        {

              fscanf(file[j], "  %s %s %s %s %s %s %s\n", &str,&str1, &str2,&str3,&str4,&str5,&str6);

             testclass1 = atof(str5);
             testclass2 = atof(str6);

            classtotal[i1][6] = round(testclass1);
            classtotal[i1][7] = round(testclass2);

            i1 =i1+1;

        }

        fclose(file[j]);
    }

    char arraytiff[50][2024];//maximum number of images is 50

   // * I AM DETECTING THE NUMBER OF MULTISPECTRAL IMAGES AS THERE COULD BE 10 , 16, 23 , ETC MULTISPECTRAL IMAGES

    DIR * dp2;
    struct dirent* ep2;
    char* text1;
    char pattern1[6];
    char* text22;
    pattern1[0]='\0';
    pattern1[1]='\0';
    pattern1[2]='\0';
    pattern1[3]='\0';
    pattern1[4]='\0';
    pattern1[5]='\0';
    text1=".tif";
    dp2 = opendir ("./");
    int acindex1, sn1;
    acindex1 = 0;

    while (ep2 = readdir (dp2))
    {
        text22 = ep2->d_name;

        sn1 = strlen(text22);;//sizeof(text22);

        pattern1[0] = text22[sn1-4];//ep->d_name[0];

        pattern1[1] = text22[sn1-3];//ep->d_name[1];

        pattern1[2] = text22[sn1-2];//ep->d_name[2];

        pattern1[3] = text22[sn1-1];//ep->d_name[3];

        pattern1[4] ='\0';
        if( strcmp(text1, pattern1) == 0  )
        {

            strcpy(arraytiff[acindex1],ep2->d_name);
            printf("%s\n",arraytiff[acindex1]);
            acindex1 = acindex1 +1;

        }
    }

    char tiffext[4];

    char namefile[40];

    strncpy(namefile,arraytiff[0],strlen(arraytiff[0])-4);

     printf(" TEST 3333333 \n ");

    // *  I AM COPYING THE PATH TO THE MULTISPECTRAL IMAGES

    (void) closedir (dp2);

    char cwd1[1024],cwdtemp1[1024],*dfg;
    getcwd(cwd1, sizeof(cwd1));

    strcat(cwd1,"/");
    char arrayimage[acindex1][2024];

    char *arrayimage2[acindex1];

    printf("%s\n",cwd);
    printf("%s\n",cwd1);

    for (int j=0; j < acindex1; j=j+1)
    {

        strcpy(cwdtemp1,cwd1);

        strcat(cwdtemp1,arraytiff[j]);

        strcpy(arrayimage[j],cwdtemp1);

        printf("%s \n",arrayimage[j]);

    }

    double data_matrix[sumtotalnrclass][acindex1];

    //#pragma omp parallel for private(j,i)
    for (int i=0; i < sumtotalnrclass; i=i+1)
    {
        for (int j=0; j < acindex1; j=j+1)
        {
            data_matrix[i][j] = 0;
        }
    }

   std::vector<Mat> image(acindex1);

    std::vector<Mat> various_images;

   // * I AM READING THE IMAGES AND I PUT THEM IN THE VECTOR OF IMAGES CALLED VARIOUS_IMAGES

    for (int i=0; i < acindex1; i=i+1)
    {

            image[i] = imread(arrayimage[i], 0);

            various_images.push_back(image[i]);

    }

    s1 = image[0].rows;
    s2 = image[0].cols;

    //int s1s2;

    s1s2 = s1*s2;

    if(! image[0].data )                              // Check for invalid input
    {
        cout <<  "Could not open or find the image" << std::endl ;
        // return -1;
    }

    // p o u b

    int *grouping_vector = new int[sumtotalnrclass] ;

    for (int k = 0; k < sumtotalnrclass; k++) {
        grouping_vector[k] = 0;

    }

    int *grp = new int[sumtotalnrclass] ;

    for (int k = 0; k < sumtotalnrclass; k++) {
        grp[k] = 0;

    }

    // * FROM THE MULTISPECTRAL IMAGES I AM SELECTING THE POINTS BASED ON THE CLASS FILES

    int indexonc;
    for (int j = 0; j < acindex1; j++) {

        indexonc = 0;
        for (int k = 0; k < acindex; k++) {

            for (int ss = 0; ss < nrclass[k]; ss++)
            {

                data_matrix[indexonc][j] =  (double) various_images[j].at<uchar>(classtotal[indexonc][7]-1,classtotal[indexonc][6]-1);

                indexonc = indexonc +1 ;

            }
        }
    }

       // * I AM COUNTING THE POINTS FROM EACH CLASS

    indexonc = 0;
    for (int k = 0; k < acindex; k=k+1) {

        for (int ss = 0; ss < nrclass[k]; ss=ss+1)
        {
            grouping_vector[indexonc] = 1+k;
            grp[indexonc]=1+k;
            indexonc = indexonc +1 ;

        }
    }

    printf("  \n ");

    double X[sumtotalnrclass][acindex1];

    for (int i = 0; i < sumtotalnrclass; i++) {
        for (int j = 0; j < acindex1; j++) {

            X[i][j] = data_matrix[i][j];

        }

    }

    int N,D,K;

    float nk[acindex1];
    float alpha,xmg[acindex1],xmk[acindex][acindex1];
    int index;
    float xdiff[sumtotalnrclass], transxdiff[sumtotalnrclass][1];
    float sum[acindex1][acindex1], sumare;

    // * I AM IMPLEMENTING THE CANONICAL VARIATES ANALYSIS METHOD BASED ON THE MATLAB VERSION OF THIS METHOD

    N = sumtotalnrclass;//sizeof(X);
    D = acindex1;//sizeof(X[0][0])-1;

    for (j=0; j < acindex1; j=j+1)
    {
        xmg[j] = 0;

    }

    K = grp[0];
    for (int j=1;j<sumtotalnrclass;j++)
    {
        if (K < grp[j])
        {
            K = grp[j];

        }
            }

    index=0;

      for (j=0;j<acindex1;j++)
    {
        alpha =0;
        for (i=0;i<sumtotalnrclass;i++)
        {
            alpha = alpha + X[i][j];

        }
        xmg[j] = alpha/sumtotalnrclass;

    }

    for (i=0; i < K; i=i+1)
    {
        for (j=0; j < acindex1; j=j+1)
        {
            xmk[i][j] = 0;

        }
    }

    float **difference = new float*[acindex1];
    for (int i = 0; i < acindex1; i++)
    {
        difference[i] = new float[acindex1];
    }

    float **difference1 = new float*[acindex1];
    for (int i = 0; i < acindex1; i++)
    {
        difference1[i] = new float[acindex1];
    }

    float suma = 0;

    double Xinter[sumtotalnrclass][acindex1];

    //printf("Hello from thread %d, nthreads %d\n", omp_get_thread_num(), omp_get_num_threads());

    //#pragma omp parallel for

   /*#pragma omp parallel for
    for(int n=0; n<10; n++)
    {
        printf(" %d", n);
    }
    printf(".\n");*/

    //#pragma omp parallel num_threads(5)
    //{

    //#pragma omp parallel for

    //#pragma omp parallel for private(k,i,c)
    for (k=1;k < acindex+1; k++)
    {
        index = 0;

        for (i=0;i<sumtotalnrclass;i++)
        {
            if (grp[i] == k)
            {

                for (int c=0;c<acindex1;c++)
                {
                    Xinter[index][c] = X[i][c];

                }
                index = index +1;
                           }

        }
        //printf("Hello from thread %d, nthreads %d\n", omp_get_thread_num(), omp_get_num_threads());

        for (int c=0;c<acindex1;c++)
        {
            suma = 0;
            for (int d=0;d<nrclass[k-1];d++)
                suma = suma + Xinter[d][c];
            xmk[k-1][c] = suma/nrclass[k-1];

        }

    }

    //}

    float **S = new float*[D];
    for (int i = 0; i < D; i++)
    {
        S[i] = new float[D];
    }

    float **B = new float*[D];
    for (int i = 0; i < D; i++)
    {
        B[i] = new float[D];
    }

    float **xsumdiff = new float*[D];
    for (int i = 0; i < D; i++)
    {
        xsumdiff[i] = new float[D];
    }

    for (i=0; i < D; i=i+1)
    {
        for (j=0; j < D; j=j+1)
        {
            S[i][j] = 0;
        }
    }

    for (i=0; i < D; i=i+1)
    {
        for (j=0; j < D; j=j+1)
        {
            xsumdiff[i][j] = 0;
        }
    }

    //#pragma omp parallel for private(i,j,c)
    for (int i=0;i<sumtotalnrclass;i++)
    {

        for (int j=0;j<acindex1;j++)
        {

            index = grp[i];
            xdiff[j] =  X[i][j] - xmk[index-1][j] ;

        }

        for (int c = 0; c < acindex1; c++)
        {
            transxdiff[c][1] = xdiff[c];

        }

        //multiply xdiff'*xdiff
        for (int c=0;c<acindex1;c++)
        {
            for (int d = 0; d < acindex1; d++)
            {
                xsumdiff[c][d] = transxdiff[c][1]*xdiff[d];
                //printf ("%f \n", xsumdiff[c][d]);
            }
        }

        for (int c = 0; c < acindex1; c++)
        {
            for (int d = 0; d < acindex1; d++)
            {
                //printf ("%f \n", S[c][d]);
                S[c][d] = S[c][d] + xsumdiff[c][d];
            }
        }

    }

    //#pragma omp parallel for private(i,j)
    for (i=0; i < acindex1; i=i+1)
    {

        for (j=0; j < acindex1; j=j+1)
        {

            S[i][j] = S[i][j]/(sumtotalnrclass-acindex);

        }
      }

    printf("\n");

    //#pragma omp parallel for private(i,j)
    for (int i=0; i < D; i=i+1)
    {
        for (int j=0; j < D; j=j+1)
        {

            B[i][j] = 0;
            //printf("%f \n",B[i][j]);

        }
    }

    double Xtemp[sumtotalnrclass][acindex1];

    //#pragma omp parallel for private(i)
    for(i = 0; i < sumtotalnrclass; i++)
    {
        //    //printf("%3d: ", i);
        for(j = 0; j < acindex1; j++)
            Xtemp[i][j] = X[i][j];

    }

    //nrclass  acindex

    //#pragma omp parallel for private(k,i,d,c)
    for (k =0;k<acindex;k++)
    {
        for (i=0;i<acindex1;i++)

        {

            difference[1][i] = nrclass[k]*(xmk[k][i] - xmg[i]);
            difference1[1][i] = xmk[k][i] - xmg[i];

        }
        for(int d = 0 ; d < acindex1 ; d++ )
            transxdiff[d][1] = difference[1][d];

        for (int d = 0; d < acindex1; d++)
        {
            for (int c = 0; c < acindex1; c++)
            {

                sum[c][d] = transxdiff[c][1]*difference1[1][d];
                //printf("%f \n",sum[c][d]);

            }

        }

        for (int c=0; c < acindex1; c=c+1)
        {
            for (int d=0; d < acindex1; d=d+1)
            {
                B[c][d] = B[c][d]+sum[c][d];

                //
                printf("%f ",B[c][d]);
            }
             printf(" \n");

        }

        printf(" \n");

    }

    double coefficient;
    double coefficient1;

    coefficient1 = acindex-1;

    coefficient = acindex/coefficient1;

     //printf(" %d \n",acindex);
    //printf(" %lf \n",coefficient);
    //printf(" %d \n",sumtotalnrclass);

    //#pragma omp parallel for private(d,c)
    for (int c=0; c < acindex1; c=c+1)
    {
        for (int d=0; d < acindex1; d=d+1)
        {

            B[c][d] = coefficient*B[c][d]/sumtotalnrclass;

            //printf("%f ",B[c][d]);

        }
        //printf(" \n");
    }

        float vectork[acindex1];

    float finalk =0;

    gsl_matrix * AA , *SS;

    gsl_vector * eval;
    gsl_matrix * evec;
    gsl_eigen_gensymmv_workspace  * w;
    //, gsl_vector_complex * eval, gsl_matrix_complex * evec, gsl_eigen_nonsymmv_workspace * w
    int err;
    err = 1;

    // * I AM USING THE GSL GNU FUNCTIONS GSL_MATRIX_ALLOC , GSL_EIGEN_GENSYMMV_ALLOC , ETC TO CALCULATE THE EIGEN VECTORS

    AA = gsl_matrix_alloc(acindex1, acindex1);
    SS = gsl_matrix_alloc(acindex1, acindex1);
    eval = gsl_vector_alloc(acindex1);
    evec  = gsl_matrix_alloc(acindex1, acindex1);
    w = gsl_eigen_gensymmv_alloc(acindex1);

    //#pragma omp parallel for private(i,j)
    for (int i = 0; i < acindex1; i++)
        for (int j = 0; j < acindex1; j++)
            gsl_matrix_set (AA, i, j, B[i][j]);// B[i][j]);

    for (int i = 0; i < acindex1; i++)
        for (int j = 0; j < acindex1; j++)
            gsl_matrix_set (SS, i, j, S[i][j]);// B[i][j]);

    int status;

    gsl_set_error_handler_off();

    status = gsl_eigen_gensymmv( AA,  SS,  eval, evec,  w) ;

    double a;
    int indexa[acindex1];

    for (i = 0; i < acindex1; i++)
    {
        indexa[i] = i;
    }

    for (i = 0; i < acindex1; i++)
    {
        indexa[i]= i;
    }

    /*for(i = 0; i < 23; i++)
     {
     printf("%3d %f %i \n", i,  gsl_vector_get(eval, i) ,indexa[i]);
     }*/

    int prst;
    //#pragma omp parallel for private(i,j)
    for (i = 0; i < acindex1; i++)
    {
        //indexa[i]= i;
        for (j = i ; j < acindex1; j++)
        {
            if (gsl_vector_get(eval, i) < gsl_vector_get(eval, j))
            {

                a = gsl_vector_get(eval, i);

                gsl_vector_set(eval, i, gsl_vector_get(eval, j) );
                prst = indexa[i];
                indexa[i]= indexa[j];
                indexa[j]= prst;

                gsl_vector_set(eval, j, a);

            }
            else
            {
                //gsl_vector_set(eval, i, (-1)*gsl_vector_get(eval, i) );
            }
        }
    }

      printf("  \n"  );
    printf(" EIGENVALUES \n ");

    for(i = 0; i < acindex1; i++)
    {
        printf("%3d %f %i \n", i,  gsl_vector_get(eval, i) ,indexa[i]);
    }

    double summ = 0;
    double Xresult[sumtotalnrclass][acindex1];

    double tempcoef[acindex1][acindex1];

    for (i = 0; i < acindex1; i++)
    {
        for (j = 0; j < acindex1; j++)
        {
            tempcoef[i][j] = gsl_matrix_get(evec,i,j) ;

        }
    }
       printf("  \n"  );
     printf(" EIGEN VECTORS \n ");

    for(i = 0; i < acindex1; i++)
    {
        //printf("%3d: ", i);
        for(j = 0; j < acindex1; j++)
            printf("%lf ",   gsl_matrix_get(evec,i,j)  );
        printf("\n");
    }

    double coef[acindex1][acindex1];

    for(i = 0; i < acindex1; i++)
    {
        //printf("%3d: ", i);
        for(j = 0; j < acindex1; j++)
            coef[i][j] =  gsl_matrix_get(evec,i,j) ;

    }

    // *  I AM CALCULATING THE NEW IMAGE POINTS BASED ON THE PARAMETERS CALCULATED BY THE GSL GNU FUNCTIONS

    //#pragma omp parallel for
    //#pragma omp parallel for private(k,i,j)
    for(i=0;i<sumtotalnrclass;i++)
    {
        for(j=0;j<acindex1;j++)
        {
            summ=0;
            for(k=0;k<acindex1;k++)
            {

                summ=summ+Xtemp[i][k]*gsl_matrix_get(evec, k, j);

            }
            Xresult[i][j]=summ;
        }

    }

     double** newmatrix1= new double*[s1];

    for (int i = 0; i < s1; i++)
    {
        newmatrix1[i] = new double[s2];
    }

     //static double newmatrix2[5412][7216];
    double** newmatrix2= new double*[s1];
    for (int i = 0; i < s1; i++)
    {
        newmatrix2[i] = new double[s2];
    }
     //static double newmatrix3[5412][7216];

    double** newmatrix3= new double*[s1];
    for (int i = 0; i < s1; i++)
    {
        newmatrix3[i] = new double[s2];
    }

    for(i = 0; i < s1; i++)
    {
        for(j = 0; j < s2; j++)
        {
            newmatrix1[i][j] = 0;
        }
    }

    for(i = 0; i < s1; i++)
    {
        for(j = 0; j < s2; j++)
        {
            newmatrix2[i][j] = 0;
        }
    }

    for(i = 0; i < s1; i++)
    {
        for(j = 0; j < s2; j++)
        {
            newmatrix3[i][j] = 0;
        }
    }

    double coef1[acindex1][acindex1];
    double score1[sumtotalnrclass][acindex1];

    //#pragma omp parallel for private(i,j)
    for (int i=0; i < sumtotalnrclass; i=i+1)
    {
        for (int j=0; j < acindex1; j=j+1)
        {
            score1[i][j] = 0;
        }
    }

    for (int i=0; i < acindex1; i=i+1)
    {
        for (int j=0; j < acindex1; j=j+1)
        {
            coef1[i][j] = 0;
        }
    }

    for (int i=0; i < sumtotalnrclass; i=i+1)
    {
        for (int j=0; j < acindex1; j=j+1)
        {
            score1[i][j] = Xresult[i][j];
        }
    }

    for (int i=0; i < acindex1; i=i+1)
    {
        for (int j=0; j < acindex1; j=j+1)
        {
            coef1[i][j] = coef[i][j];
        }
    }

     // *  I AM IMPLEMENTING THE POST-PROCESSING WHICH EXIST IN THE ORIGINAL CANONICAL VARIATES ANALYSIS METHOD IMPLEMENTED INITIALLY IN MATLAB

    //#pragma omp parallel for private(j,i,k)
    for(i = 0; i < s1; i++)
    {
        for(j = 0; j < s2; j++)
        {
            for(k = 0; k < acindex1; k++)
            {
                newmatrix1[i][j] = newmatrix1[i][j]  + ( (double)various_images[k].at<uchar>(i,j) )*coef1[k][0];

                newmatrix2[i][j] = newmatrix2[i][j]  + ( (double)various_images[k].at<uchar>(i,j) )*coef1[k][1];

                newmatrix3[i][j] = newmatrix3[i][j]  + ( (double)various_images[k].at<uchar>(i,j) )*coef1[k][2];

            }

        }

    }

    double minimum1,minimum2,minimum3;
    double maximum1,maximum2,maximum3;

    double **range_cv1 = new double*[s1];
    for (int i = 0; i < s1; i++)
    {
        range_cv1[i] = new double[s2];
    }

    double **range_cv2 = new double*[s1];
    for (int i = 0; i < s1; i++)
    {
        range_cv2[i] = new double[s2];
    }

    double **range_cv3 = new double*[s1];
    for (int i = 0; i < s1; i++)
    {
        range_cv3[i] = new double[s2];
    }

    minimum1= newmatrix1[0][0];
    for(int i=0; i<s1; i++)
    {
        for(int j=0; j<s2; j++)
        {
            if( newmatrix1[i][j]< minimum1)
            {
                minimum1=newmatrix1[i][j];
            }
        }
    }

    minimum2= newmatrix2[0][0];
    for(int i=0; i<s1; i++)
    {
        for(int j=0; j<s2; j++)
        {
            if( newmatrix2[i][j]< minimum2)
            {
                minimum2=newmatrix2[i][j];
            }
        }
    }

    minimum3= newmatrix3[0][0];
    for(int i=0; i<s1; i++)
    {
        for(int j=0; j<s2; j++)
        {
            if( newmatrix3[i][j]< minimum3)
            {
                minimum3=newmatrix3[i][j];
            }
        }
    }

    maximum1= newmatrix1[0][0];
    for(int i=0; i<s1; i++)
    {
        for(int j=0; j<s2; j++)
        {
            if( newmatrix1[i][j]> maximum1)
            {
                maximum1=newmatrix1[i][j];
            }
        }
    }

    maximum2= newmatrix2[0][0];
    for(int i=0; i<s1; i++)
    {
        for(int j=0; j<s2; j++)
        {
            if( newmatrix2[i][j]> maximum2)
            {
                maximum2=newmatrix2[i][j];
            }
        }
    }

    maximum3= newmatrix3[0][0];
    for(int i=0; i<s1; i++)
    {
        for(int j=0; j<s2; j++)
        {
            if( newmatrix3[i][j]> maximum3)
            {
                maximum3=newmatrix3[i][j];
            }
        }
    }

    range_cv1 = range_map(newmatrix1, minimum1, maximum1);

    range_cv2 = range_map(newmatrix2, minimum2, maximum2);

    range_cv3 = range_map(newmatrix3, minimum3, maximum3);

    int aa,bb,cc;

     printf("  \n"  );

    printf(" A SET OF 8 PICTURES ARE PRODUCED ");

    cv::Mat img16(s1, s2, CV_8UC3);

    vector<int> compression_params;
    compression_params.push_back(CV_IMWRITE_PNG_COMPRESSION);
    compression_params.push_back(9);

    for (int i = 0; i < s1; i++) {
        for (int j = 0; j < s2; j++) {

            img16.at<cv::Vec3b>(i,j)[1] =  (255-range_cv2[i][j]);
            img16.at<cv::Vec3b>(i,j)[2] = (255-range_cv1[i][j]);
            img16.at<cv::Vec3b>(i,j)[3] = (255-range_cv3[i][j]);

        }
    }

    char fileword[100];
    strcpy(fileword,namefile);
    strcat(fileword,"_image_rangecvasecond.tiff");

    //end = clock();
    //cpu_time_used = ((double) (end - start)) / CLOCKS_PER_SEC;
    //printf("  \n"  );

    //printf(" %lf ", cpu_time_used  );

    //exit(1);

    cv::imwrite(fileword, img16);

    cv::Mat img162(s1, s2, CV_8UC3);

    for (int i = 0; i < s1; i++) {
        for (int j = 0; j < s2; j++) {

            img162.at<cv::Vec3b>(i,j)[1] =  range_cv2[i][j];
            img162.at<cv::Vec3b>(i,j)[2] = range_cv1[i][j];
            img162.at<cv::Vec3b>(i,j)[3] = range_cv3[i][j];

        }
    }

    char fileword1[100];
    strcpy(fileword1,namefile);
    strcat(fileword1,"_image_rangecva.tiff");

     cv::imwrite(fileword1, img162);

    double minscore1,minscore2,minscore3;
    double maxscore1,maxscore2,maxscore3;

    minscore1=0;
    minscore1= score1[0][0];
    for(int i=0; i<sumtotalnrclass; i++)
    {

        if( score1[i][0]< minscore1)
        {
            minscore1=score1[i][0];
        }

    }

    minscore2==0;
    minscore2= score1[0][1];
    for(int i=0; i<sumtotalnrclass; i++)
    {

        if( score1[i][1]< minscore2)
        {
            minscore2=score1[i][1];
        }

    }

    minscore3==0;
    minscore3= score1[0][2];
    for(int i=0; i<sumtotalnrclass; i++)
    {

        if( score1[i][2]< minscore3)
        {
            minscore3=score1[i][2];
        }

    }

    maxscore1=0;
    maxscore1= score1[0][0];
    for(int i=0; i<sumtotalnrclass; i++)
    {

        if( score1[i][0]> maxscore1)
        {
            maxscore1=score1[i][0];
        }

    }

    maxscore2=0;
    maxscore2= score1[0][1];
    for(int i=0; i<sumtotalnrclass; i++)
    {

        if( score1[i][1]> maxscore2)
        {
            maxscore2=score1[i][1];
        }

    }

    maxscore3=0;
    maxscore3= score1[0][2];
    for(int i=0; i<sumtotalnrclass; i++)
    {

        if( newmatrix3[i][2]> maxscore3)
        {
            maxscore3=score1[i][2];
        }

    }

    double **score_cv1 = new double*[s1];
    for (int i = 0; i < s1; i++)
    {
        score_cv1[i] = new double[s2];
    }

    double **score_cv2 = new double*[s1];
    for (int i = 0; i < s1; i++)
    {
        score_cv2[i] = new double[s2];
    }

    double **score_cv3 = new double*[s1];
    for (int i = 0; i < s1; i++)
    {
        score_cv3[i] = new double[s2];
    }

    score_cv1 = range_map(newmatrix1, minscore1, maxscore1);

    score_cv2 = range_map(newmatrix2, minscore2, maxscore2);

    score_cv3 = range_map(newmatrix3, minscore3, maxscore3);

    cv::Mat img163(s1, s2, CV_8UC3);

    for (int i = 0; i < s1; i++) {
        for (int j = 0; j < s2; j++) {

            img163.at<cv::Vec3b>(i,j)[1] =  score_cv2[i][j];
            img163.at<cv::Vec3b>(i,j)[2] = score_cv1[i][j];
            img163.at<cv::Vec3b>(i,j)[3] = score_cv3[i][j];

        }
    }

    char fileword2[100];
    strcpy(fileword2,namefile);
    strcat(fileword2,"_image_score_rangecva.tiff");

    cv::imwrite(fileword2, img163);

    cv::Mat img164(s1, s2, CV_8UC3);

     for (int i = 0; i < s1; i++) {
        for (int j = 0; j < s2; j++) {

            img164.at<cv::Vec3b>(i,j)[1] =  (255-score_cv2[i][j]);
            img164.at<cv::Vec3b>(i,j)[2] = (255-score_cv1[i][j]);
            img164.at<cv::Vec3b>(i,j)[3] = (255-score_cv3[i][j]);

        }
    }

    char fileword3[100];
    strcpy(fileword3,namefile);
    strcat(fileword3,"_image_score_rangecvasecond.tiff");

    cv::imwrite(fileword3, img164);

       // * I AM PRODUCING THE 8 BIT IMAGES OF THE PERCENTILE

    double required_percentiles[8];

    required_percentiles[0] = 0.01;
    required_percentiles[1] = 0.1;
    required_percentiles[2] = 1;
    required_percentiles[3] = 5;
    required_percentiles[4] = 99.99;
    required_percentiles[5] = 99.9;
    required_percentiles[6] = 99;
    required_percentiles[7] = 95;

    int sizerows,sizecolumns;

    int n_percentiles;

    n_percentiles = (int) 8 / 2;

    double *percentiles_cv1 = new double[5];
    double *percentiles_cv2 = new double[5];
    double *percentiles_cv3 = new double[5];
    double *percentiles_cv4 = new double[5];
    double *percentiles_cv5 = new double[5];

    for (int i = 0; i < 5; i++)
    {
        percentiles_cv1[i] = 0;
        percentiles_cv2[i] = 0;
        percentiles_cv3[i] = 0;
        percentiles_cv4[i] = 0;
        percentiles_cv5[i] = 0;

    }

    percentiles_cv1 = Percentile(newmatrix1, required_percentiles);
    percentiles_cv2 = Percentile(newmatrix2, required_percentiles);
    percentiles_cv3 = Percentile(newmatrix3, required_percentiles);

    double **per_cv1 = new double*[s1];
    for (int i = 0; i < s1; i++)
    {
        per_cv1[i] = new double[s2];
    }

    double **per_cv2 = new double*[s1];
    for (int i = 0; i < s1; i++)
    {
        per_cv2[i] = new double[s2];
    }

    double **per_cv3 = new double*[s1];
    for (int i = 0; i < s1; i++)
    {
        per_cv3[i] = new double[s2];
    }

    cv::Mat img165(s1, s2, CV_8UC3);
    char str22[35];

    for (int i = 0; i < n_percentiles ; i++)
    {
        per_cv1 = range_map(newmatrix1,percentiles_cv1[i] , percentiles_cv1[i+n_percentiles]);

        per_cv2 = range_map(newmatrix2, percentiles_cv2[i], percentiles_cv2[i+n_percentiles]);

        per_cv3 = range_map(newmatrix3, percentiles_cv3[i], percentiles_cv3[i+n_percentiles]);

        for (int i = 0; i < s1; i++) {
            for (int j = 0; j < s2; j++) {

                img165.at<cv::Vec3b>(i,j)[1] =  per_cv2[i][j];
                img165.at<cv::Vec3b>(i,j)[2] =  per_cv1[i][j];
                img165.at<cv::Vec3b>(i,j)[3] =  per_cv3[i][j];

            }
        }

        sprintf(str22,"_image_percentile_%lfcva.tiff", required_percentiles[i]);

        strcpy(fileword3,namefile);
        strcat(fileword3,str22);

        cv::imwrite(fileword3, img165);

    }

    end = clock();
    cpu_time_used = ((double) (end - start)) / CLOCKS_PER_SEC;

    printf(" %lf ", cpu_time_used  );

    return;
}

//double* Percentile(double arr[5412][7216], double percentile[8])
double* Percentile(double** arr, double percentile[8])
{

    int index,size1;

    index =0;

    double *resultingvector = new double[s1s2];

    size1 = sizeof(percentile);

    for(int i = 0; i < s1; i++)
    {
        for(int j = 0; j < s2; j++)
        {

            resultingvector[index] = arr[i][j];
            index = index +1;
        }

    }

    double aa;

    gsl_vector * v = gsl_vector_alloc (s1s2);

    for (int i = 0; i < s1s2; i++)
    {
        gsl_vector_set (v, i,resultingvector[i] );
    }

    gsl_sort_vector(v);

    for (int i = 0; i < s1s2; i++)
    {
        resultingvector[i] = gsl_vector_get (v, i);
    }

    int index2;

    double *percentile_values = new double[size1];

    for(int i = 0; i < size1; i++)
    {
        index2 = round(1 + (percentile[i]/100)*(s1s2-1));
        percentile_values[i] = resultingvector[index2];
    }

    return percentile_values;
}

int compare_doubles (const double * a,
                                 const double * b)
{
    if (*a > *b)
        return 1;
    else if (*a < *b)
        return -1;
    else
        return 0;
}

double* quickSort( double a[], long l, long r)
{
    long j;

    //r = 39052992;
    if( l < r )
    {
        // divide and conquer
        j = partition( a, l, r);
        a = quickSort( a, l, j-1);
        a = quickSort( a, j+1, r);
    }
    //for(int i = 0; i < 39052992; i++)
    //{
    printf(" %d  \n",  j  );
    //}
    return a;
}

long partition(double a[], long l, long r) {
    long  i, j;
    double pivot,t;
    pivot = a[l];
    i = l; j = r+1;

    while( 1)
    {
        do ++i; while( a[i] <= pivot && i <= r );
        do --j; while( a[j] > pivot );
        if( i >= j ) break;
        t = a[i]; a[i] = a[j]; a[j] = t;
    }
    t = a[l]; a[l] = a[j]; a[j] = t;
    return j;
}

//double** range_map(double newmatrix[5412][7216], double minimum, double maximum)
double** range_map(double** newmatrix, double minimum, double maximum)

{

    //static double newimage[5412][7216];

    double **newimage = new double*[s1];
    for (int i = 0; i < s1; i++)
    {
        newimage[i] = new double[s2];
    }

    double out_low, out_high, out_range_low, out_range_high,v,v2;

    out_low = 0;
    out_high =255;
    out_range_low = 0;
    out_range_high = 255;

    double in_low,in_high;

    in_low = minimum;
    in_high = maximum;

    for(int i = 0; i < s1; i++)
    {
        for(int j = 0; j < s2; j++)
        {
            v =  newmatrix[i][j];
            if (v < in_low)
                v2 = out_range_low;
            else if (v > in_high)
                v2 = out_range_high;
            else
                v2 = ((v - in_low) / (in_high - in_low)) * (out_high - out_low) + out_low;

            //newimage[i][j] = uint8(floor(v2 + 0.5));
            newimage[i][j] = uint8_t(floor(v2+0.5));

        }
    }
    return newimage;
}


text:
/*
 * Copyright 2008-2016 Jan Gasthaus
 *
 * Licensed under the Apache License, Version 2.0 (the "License");
 * you may not use this file except in compliance with the License.
 * You may obtain a copy of the License at
 *
 *     http://www.apache.org/licenses/LICENSE-2.0
 *
 * Unless required by applicable law or agreed to in writing, software
 * distributed under the License is distributed on an "AS IS" BASIS,
 * WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied.
 * See the License for the specific language governing permissions and
 * limitations under the License.
 */

#ifndef UTILS_H_
#define UTILS_H_

#include <string>
#include <cmath>
#include <map>
#include <iostream>
#include <fstream>
#include <sstream>
#include <cassert>
#include <algorithm>
#include <gsl/gsl_sf_gamma.h>

////////////////////////////////////////////////////////////////////////////////
/////////////////   SOME USEFUL MACROS   ///////////////////////////////////////
////////////////////////////////////////////////////////////////////////////////
#ifdef DEBUG
#define tracer if (0) ; else std::cerr
#define DBG if (1)
#else
#define tracer if (1) ; else std::cerr
#define DBG if (1) ; else
#endif

  // some syntactic sugar for implementing interfaces using (multiple) inheritance
#define interface class
#define implements public

  // A macro to disallow the copy constructor and operator= functions
  // This should be used in the private: declarations for a class
#define DISALLOW_COPY_AND_ASSIGN(TypeName) \
  TypeName(const TypeName&);               \
  void operator=(const TypeName&)

namespace gatsby { namespace libplump {

/**
 * Alias: std::vector<double>
 */
typedef std::vector<double> d_vec;

/**
 * Alias: std::vector<std::vector<double> >
 */
typedef std::vector<d_vec> d_vec_vec;

/**
 * Alias: std::vector<unsigned int>
 */
typedef std::vector<unsigned int> ui_vec;

/**
 * Alias: std::vector<std::vector<unsigned int> >
 */
typedef std::vector<ui_vec> ui_vec_vec;

/**
 * Compute log2(x) -- the log2 provided by GCC is somewhat weird.
 */
inline double log2(double x){
  static const double LOG2 = std::log(2.);
  return std::log(x)/LOG2;
}

/**
 * Read items of type E from a file and push them into a sequence of type
 * S (that has to support push_back(E item).
 *
 * Items are read from the file sizeof(E) bytes at a time. No error checking
 * is performed.
 *
 * @tparam E Type of element to read from file, elements will be read sizeof(E)
 *           bytes at a time.
 * @tparam S Type of sequence to push items into;
 *           must support push_back(E item)
 * @param fileName Name of the file to read from
 * @param seq Sequence to push items into
 * @param limit Maximum number of items to push (default: 0, no limit)
 */
template<typename E, class S>
  inline void pushFileToVec(const std::string& fileName, S& seq, int limit = 0){
    std::ifstream in;
    in.open(fileName.c_str(),std::ios::binary | std::ios_base::in);
    char buffer[sizeof(E)];
    in.read(buffer,sizeof(E));
    int j = 0;
    if (limit==0)
      limit = -1;
    while (!in.eof() && j!=limit){
      seq.push_back(*((E*)buffer));
      in.read(buffer,sizeof(E));
      j++;
    }
    in.close();
  }

/**
 * Push each character of string s individually into container seq of type S,
 * which has to support push_back(char c);
 *
 * @tparam S type of container to push into
 * @param s input string to push
 * @param seq output container to push into
 */
template<class S>
  inline void pushStringToVec(const std::string& s, S& seq){
    for (size_t i=0;i<s.length();i++) {
      seq.push_back(s[i]);
    }
  }

/**
 * Convert a map<T,T> where T is an integer type with
 * range 0...MAX to a vector v of length MAX+1 so that
 * v[x] = m[x] for x=1,...,MAX. v[x] = 0 if m.count(x) = 0.
 *
 * @tparam T integer type
 * @param m input map
 * @param v output vector
 * @returns pointer to histogram vector on the heap
 */
template<typename T>
  inline void mapHistogramToVectorHistogram(const std::map<T,T>& m, std::vector<T>& v) {
    T max = (*max_element(m.begin(), m.end())).first;
    v.clear();
    v.assign(max + 1,0);
    for (typename std::map<T,T>::iterator i = m.begin(); i != m.end(); ++i) {
      v[(*i).first] = (*i).second;
    }
    return v;
  }

/**
 * Count the number of times the elements k=0...max-1 occur in vec and
 * return the result in ret[k].
 */
template<typename T>
  inline std::vector<T> vec2hist(std::vector<T>& vec, T max) {
    std::vector<T> ret(max,0);
    for (unsigned int i = 0; i < vec.size(); ++i) {
      ++ret[vec[i]];
    }
    return ret;
  }

/**
 * Compute the mean of the elements in a sequence.
 *
 * @tparam T iteratable sequence type; must have a const_iterator member
 *           as well as begin() and end() methods.
 * @param in sequence to compute the mean of
 * @returns the mean of the elements in the sequence
 */
template<typename T>
  inline double mean(const T& in) {
    double mean = 0.0;
    size_t j = 0;
    for (typename T::const_iterator i = in.begin(); i != in.end(); ++i) {
      mean += (((double)*i)-mean)/(++j);
    }
    return mean;
  }

/**
 * Sum the elements in a sequence.
 *
 * @tparam T iteratable sequence type; must have a const_iterator member
 *           as well as begin() and end() methods.
 * @param in sequence to compute the sum of
 * @returns the sum of the elements in the sequence
 */
template<typename T>
  inline typename T::value_type sum(const T& in) {
    typename T::value_type sum = 0;
    for (typename T::const_iterator i = in.begin(); i != in.end(); ++i) {
      sum += *i;
    }
    return sum;
  }

inline void log2_vec(d_vec& in) {
  for (size_t i=0;i<in.size();i++){
    in[i] = log2(in[i]);
  }
}

inline void exp_vec(d_vec& in) {
  for (size_t i=0;i<in.size();i++){
    in[i] = exp(in[i]);
  }
}

/**
 * Multiply all elements of a vector by a constant.
 */
template<typename elem_t>
  inline void mult_vec(std::vector<elem_t>& in, elem_t mult) {
    for (size_t i=0;i<in.size();i++){
      in[i] = in[i] * mult;
    }
  }

/**
 * Add a constant to all elements of a vector.
 */
template<typename elem_t>
  inline void add_vec(std::vector<elem_t>& in, elem_t add) {
    for (size_t i=0;i<in.size();i++){
      in[i] = in[i] + add;
    }
  }

/**
 * Subtract max element fromall elements of a vector.
 */
template<typename elem_t>
  inline void subMax_vec(std::vector<elem_t>& in) {
    elem_t max = *std::max_element(in.begin(), in.end());
    // if (max == -INFINITY) {
    //   for (size_t i=0;i<in.size();i++){
    //     in[i] = 0;
    //   }
    // }
    for (size_t i=0;i<in.size();i++){
      in[i] = in[i] - max;
    }
  }

/**
 * Add two vectors elementwise.
 */
template<typename elem_t>
  inline void add_vec(std::vector<elem_t>& inout, const std::vector<elem_t>& add) {
    std::transform(inout.begin(), inout.end(), add.begin(), inout.begin(),
                  std::plus<elem_t>());
  }

/**
 * Elementwise multiplication. Results is returned in the first argument.
 */
template<typename elem_t>
  inline void mult_vec(std::vector<elem_t>& inout, std::vector<elem_t> in) {
    assert(inout.size() == in.size());
    for (size_t i=0;i<inout.size();i++){
      inout[i] *= in[i];
    }
  }

/**
 * Average the values in multiple vectors of the same length.
 **/
inline d_vec average(const d_vec_vec& ins) {
  d_vec out(ins[0].size(), 0);
  for (int j = 0; j < ins.size(); j++) {
    add_vec(out, ins[j]);
  }
  mult_vec(out, 1./ins.size());
  return out;
}

template<typename elem_t>
  inline double prob2loss(std::vector<elem_t> in) {
    double out = 0.0;
    for (size_t i=0;i<in.size();i++){
      out -= log2(in[i]);
    }
    return out/in.size();
  }

inline bool closeTo(double value, double target, double tol=10e-4) {
  return std::abs(value-target) < tol;
}

template<typename Iterable>
  inline std::string iterableToString(const Iterable input) {
    std::ostringstream output;
    for(typename Iterable::const_iterator it = input.begin(); it!=input.end();it++) {
      output << *it << ", ";
    }
    return output.str();
  }

template<typename Iterable>
  void iterableToCSVFile(const Iterable input, std::string fn) {
    std::ofstream output(fn.c_str());
    output.precision(10);
    for(typename Iterable::const_iterator it = input.begin(); it!=input.end();it++) {
      output << *it << ", ";
    }
  }

/**
 * computes log(exp(a) + exp(b)) while avoiding numerical
 * instabilities.
 */
inline double logsumexp(double a, double b) {
  // choose c to be the one that is largest in abs value
  double c = (a>b)?a:b;
  return (log(exp(a-c) + exp(b-c))) + c;
}

inline double sigmoid(double x) {
  return 1.0/(1.0 + exp(-x));
}

// logit = inverse sigmoid
inline double logit(double x) {
  return log(x) - log(1-x);
}

inline double logKramp(double base, double inc, double lim) {
  if (inc == 0)
      return lim*log(base);
  if (lim <= 0) {
    return 0;
  } else {
    return lim*log(inc) + gsl_sf_lnpoch(base/inc, lim);
  }
}

inline double kramp(double base, double inc, double lim) {
  return exp(logKramp(base, inc, lim));
}

static clock_t global_clock;

/*
 * Start timing
 */
inline void tic() {
  global_clock = clock();
}

/**
 * Return number of seconds since last tic()
 */
inline double toc() {
  return (clock() - global_clock)/(double)CLOCKS_PER_SEC;

}

inline std::string makeProgressBarString(double percentDone, int total=10) {
  std::ostringstream out;

  int numDone = (int)floor(total * percentDone);

  out << "[";
  for (int i=0; i<numDone; i++) {
    out << "=";
  }
  out << ">";
  for (int i=numDone; i<total; i++) {
    out << " ";
  }
  out << "]";
  return out.str();
}

}} // namespace gatsby::libplump

#endif /* UTILS_H_ */

ext: h

max_stars_repo_name: jgasthaus/libPLUMP

max_issues_repo_name: jgasthaus/libPLUMP

max_forks_repo_name: jgasthaus/libPLUMP


text:
/**
 * @file zheevd.c
 *
 *  PLASMA computational routines
 *  Release Date: November, 15th 2009
 *  PLASMA is a software package provided by Univ. of Tennessee,
 *  Univ. of California Berkeley and Univ. of Colorado Denver
 *
 * @version 2.6.0
 * @author Azzam Haidar
 * @author Hatem Ltaief
 * @date 2010-11-15
 * @precisions normal z -> s d c
 *
 **/
#include <lapacke.h>
#include "common.h"
/***************************************************************************//**
 *
 * @ingroup PLASMA_Complex64_t
 *
 *  PLASMA_zheevd - Computes all eigenvalues and, optionally,
 *  eigenvectors of a complex Hermitian matrix A. The matrix A is
 *  preliminary reduced to tridiagonal form using a two-stage
 *  approach:
 *  First stage: reduction to band tridiagonal form;
 *  Second stage: reduction from band to tridiagonal form.
 *
 *******************************************************************************
 *
 * @param[in] jobz
 *          Intended usage:
 *          = PlasmaNoVec: computes eigenvalues only;
 *          = PlasmaVec: computes eigenvalues and eigenvectors.
 *
 * @param[in] uplo
 *          Specifies whether the matrix A is upper triangular or
 *          lower triangular:
 *          = PlasmaUpper: Upper triangle of A is stored;
 *          = PlasmaLower: Lower triangle of A is stored.
 *
 * @param[in] N
 *          The order of the matrix A. N >= 0.
 *
 * @param[in,out] A
 *          On entry, the symmetric (or Hermitian) matrix A.
 *          If uplo = PlasmaUpper, the leading N-by-N upper triangular
 *          part of A contains the upper triangular part of the matrix
 *          A, and the strictly lower triangular part of A is not
 *          referenced.
 *          If uplo = PlasmaLower, the leading N-by-N lower triangular
 *          part of A contains the lower triangular part of the matrix
 *          A, and the strictly upper triangular part of A is not
 *          referenced.
 *          On exit, the lower triangle (if uplo = PlasmaLower) or the
 *          upper triangle (if uplo = PlasmaUpper) of A, including the
 *          diagonal, is destroyed.
 *
 * @param[in] LDA
 *          The leading dimension of the array A. LDA >= max(1,N).
 *
 * @param[out] W
 *          On exit, if info = 0, the eigenvalues.
 *
 * @param[in, out] descT
 *          On entry, descriptor as return by PLASMA_Alloc_Workspace_zheevd
 *          On exit, contains auxiliary factorization data.
 *
 * @param[out] Q
 *          On exit, if jobz = PlasmaVec and info = 0, the eigenvectors.
 *
 * @param[in] LDQ
 *          The leading dimension of the array Q. LDQ >= max(1,N).
 *
 *******************************************************************************
 *
 * @return
 *          \retval PLASMA_SUCCESS successful exit
 *          \retval <0 if -i, the i-th argument had an illegal value
 *          \retval >0 if INFO = i, the algorithm failed to converge; i
 *               off-diagonal elements of an intermediate tridiagonal
 *               form did not converge to zero.
 *
 *******************************************************************************
 *
 * @sa PLASMA_zheevd_Tile
 * @sa PLASMA_zheevd_Tile_Async
 * @sa PLASMA_cheevd
 * @sa PLASMA_dsyev
 * @sa PLASMA_ssyev
 *
 ******************************************************************************/
int PLASMA_zheevd(PLASMA_enum jobz, PLASMA_enum uplo, int N,
                 PLASMA_Complex64_t *A, int LDA,
                 double *W,
                 PLASMA_desc *descT,
                 PLASMA_Complex64_t *Q, int LDQ)
{
    int NB;
    int status;
    plasma_context_t *plasma;
    PLASMA_sequence *sequence = NULL;
    PLASMA_request request = PLASMA_REQUEST_INITIALIZER;
    PLASMA_desc descA;

    plasma = plasma_context_self();
    if (plasma == NULL) {
        plasma_error("PLASMA_zheevd", "PLASMA not initialized");
        return PLASMA_ERR_NOT_INITIALIZED;
    }

    /* Check input arguments */
    if (jobz != PlasmaNoVec && jobz != PlasmaVec) {
        plasma_error("PLASMA_zheevd", "illegal value of jobz");
        return -1;
    }
    if (uplo != PlasmaLower && uplo != PlasmaUpper) {
        plasma_error("PLASMA_zheevd", "illegal value of uplo");
        return -2;
    }
    if (N < 0) {
        plasma_error("PLASMA_zheevd", "illegal value of N");
        return -3;
    }
    if (LDA < max(1, N)) {
        plasma_error("PLASMA_zheevd", "illegal value of LDA");
        return -5;
    }
    if (LDQ < max(1, N)) {
        plasma_error("PLASMA_zheevd", "illegal value of LDQ");
        return -9;
    }

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

    /* Tune NB & IB depending on N; Set NBNB */
    status = plasma_tune(PLASMA_FUNC_ZHEEVD, N, N, 0);
    if (status != PLASMA_SUCCESS) {
        plasma_error("PLASMA_zheevd", "plasma_tune() failed");
        return status;
    }

    /* Set NT */
    NB = PLASMA_NB;

    plasma_sequence_create(plasma, &sequence);

    if ( PLASMA_TRANSLATION == PLASMA_OUTOFPLACE ) {
        plasma_zooplap2tile( descA, A, NB, NB, LDA, N, 0, 0, N, N, sequence, &request,
                             plasma_desc_mat_free(&(descA)) );
    } else {
        plasma_ziplap2tile( descA, A, NB, NB, LDA, N, 0, 0, N, N,
                            sequence, &request);
    }

    /* Call the tile interface */
    PLASMA_zheevd_Tile_Async(jobz, uplo, &descA, W, descT, Q, LDQ, sequence, &request);

    if ( PLASMA_TRANSLATION == PLASMA_OUTOFPLACE ) {
        plasma_zooptile2lap( descA, A, NB, NB, LDA, N,  sequence, &request);
        plasma_dynamic_sync();
        plasma_desc_mat_free(&descA);
    } else {
        plasma_ziptile2lap( descA, A, NB, NB, LDA, N,  sequence, &request);
        plasma_dynamic_sync();
    }

    status = sequence->status;
    plasma_sequence_destroy(plasma, sequence);
    return status;
}
/***************************************************************************//**
 *
 * @ingroup PLASMA_Complex64_t_Tile
 *
 *  PLASMA_zheevd_Tile - Computes all eigenvalues and, optionally, eigenvectors of a
 *  complex Hermitian matrix A using a two-stage approach:
 *  First stage: reduction to band tridiagonal form;
 *  Second stage: reduction from band to tridiagonal form.
 *
 *  Operates on matrices stored by tiles.
 *  All matrices are passed through descriptors.
 *  All dimensions are taken from the descriptors.
 *
 *******************************************************************************
 *
 * @param[in] jobz
 *          Intended usage:
 *          = PlasmaNoVec: computes eigenvalues only;
 *          = PlasmaVec: computes eigenvalues and eigenvectors.
 *
 * @param[in] uplo
 *          Specifies whether the matrix A is upper triangular or
 *          lower triangular:
 *          = PlasmaUpper: Upper triangle of A is stored;
 *          = PlasmaLower: Lower triangle of A is stored.
 *
 * @param[in,out] A
 *          On entry, the symmetric (or Hermitian) matrix A.
 *          If uplo = PlasmaUpper, the leading N-by-N upper triangular
 *          part of A contains the upper triangular part of the matrix
 *          A, and the strictly lower triangular part of A is not
 *          referenced.
 *          If UPLO = 'L', the leading N-by-N lower triangular part of
 *          A contains the lower triangular part of the matrix A, and
 *          the strictly upper triangular part of A is not referenced.
 *          On exit, if jobz = PlasmaVec, then if return value = 0, A
 *          contains the orthonormal eigenvectors of the matrix A.
 *          If jobz = PlasmaNoVec, then on exit the lower triangle (if
 *          uplo = PlasmaLower) or the upper triangle (if uplo =
 *          PlasmaUpper) of A, including the diagonal, is destroyed.*
 *
 * @param[out] W
 *          On exit, if info = 0, the eigenvalues.
 *
 * @param[in,out] T
 *          On entry, descriptor as return by
 *          PLASMA_Alloc_Workspace_zheevd
 *          On exit, contains auxiliary factorization data.
 *
 * @param[out] Q
 *          On exit, if jobz = PlasmaVec and info = 0, the eigenvectors.
 *
 * @param[in] LDQ
 *          The leading dimention of the eigenvectors matrix Q. LDQ >= max(1,N).
 *
 *******************************************************************************
 *
 * @return
 *          \retval PLASMA_SUCCESS successful exit
 *          \retval <0 if -i, the i-th argument had an illegal value
 *          \retval >0 if INFO = i, the algorithm failed to converge; i
 *               off-diagonal elements of an intermediate tridiagonal
 *               form did not converge to zero.
 *
 *******************************************************************************
 *
 * @sa PLASMA_zheevd_Tile
 * @sa PLASMA_zheevd_Tile_Async
 * @sa PLASMA_cheevd
 * @sa PLASMA_dsyev
 * @sa PLASMA_ssyev
 *
 ******************************************************************************/
int PLASMA_zheevd_Tile(PLASMA_enum jobz, PLASMA_enum uplo,
                      PLASMA_desc *A, double *W,
                      PLASMA_desc *T, PLASMA_Complex64_t *Q, int LDQ)
{
    plasma_context_t *plasma;
    PLASMA_sequence *sequence = NULL;
    PLASMA_request request = PLASMA_REQUEST_INITIALIZER;
    int status;

    plasma = plasma_context_self();
    if (plasma == NULL) {
        plasma_fatal_error("PLASMA_zheevd_Tile", "PLASMA not initialized");
        return PLASMA_ERR_NOT_INITIALIZED;
    }
    plasma_sequence_create(plasma, &sequence);
    PLASMA_zheevd_Tile_Async(jobz, uplo, A, W, T, Q, LDQ, sequence, &request);
    plasma_dynamic_sync();
    status = sequence->status;
    plasma_sequence_destroy(plasma, sequence);
    return status;
}

/***************************************************************************//**
 *
 * @ingroup PLASMA_Complex64_t_Tile_Async
 *
 *  PLASMA_zheevd_Tile_Async - Computes all eigenvalues and,
 *  optionally, eigenvectors of a complex Hermitian matrix A using a
 *  two-stage approach:
 *  First stage: reduction to band tridiagonal form;
 *  Second stage: reduction from band to tridiagonal form.
 *
 *  May return before the computation is finished.
 *  Allows for pipelining of operations at runtime.
 *
 *******************************************************************************
 *
 * @param[in] jobz
 *          Intended usage:
 *          = PlasmaNoVec: computes eigenvalues only;
 *          = PlasmaVec: computes eigenvalues and eigenvectors.
 *
 * @param[in] uplo
 *          Specifies whether the matrix A is upper triangular or
 *          lower triangular:
 *          = PlasmaUpper: Upper triangle of A is stored;
 *          = PlasmaLower: Lower triangle of A is stored.
 *
 * @param[in,out] A
 *          On entry, the symmetric (or Hermitian) matrix A.
 *          If uplo = PlasmaUpper, the leading N-by-N upper triangular
 *          part of A contains the upper triangular part of the matrix
 *          A, and the strictly lower triangular part of A is not
 *          referenced.
 *          If UPLO = 'L', the leading N-by-N lower triangular part of
 *          A contains the lower triangular part of the matrix A, and
 *          the strictly upper triangular part of A is not referenced.
 *          On exit, if jobz = PlasmaVec, then if return value = 0, A
 *          contains the orthonormal eigenvectors of the matrix A.
 *          If jobz = PlasmaNoVec, then on exit the lower triangle (if
 *          uplo = PlasmaLower) or the upper triangle (if uplo =
 *          PlasmaUpper) of A, including the diagonal, is destroyed.*
 *
 * @param[out] W
 *          On exit, if info = 0, the eigenvalues.
 *
 * @param[in,out] T
 *          On entry, descriptor as return by
 *          PLASMA_Alloc_Workspace_zheevd
 *          On exit, contains auxiliary factorization data.
 *
 * @param[out] Q
 *          On exit, if jobz = PlasmaVec and info = 0, the eigenvectors.
 *
 * @param[in] LDQ
 *          The leading dimention of the eigenvectors matrix Q. LDQ >= max(1,N).
 *
 * @param[in] sequence
 *          Identifies the sequence of function calls that this call belongs to
 *          (for completion checks and exception handling purposes).
 *
 * @param[out] request
 *          Identifies this function call (for exception handling purposes).
 *
 *******************************************************************************
 *
 * @sa PLASMA_zheevd
 * @sa PLASMA_zheevd_Tile
 * @sa PLASMA_cheevd_Tile_Async
 * @sa PLASMA_dsyev_Tile_Async
 * @sa PLASMA_ssyev_Tile_Async
 *
 ******************************************************************************/
int PLASMA_zheevd_Tile_Async(PLASMA_enum jobz, PLASMA_enum uplo,
                            PLASMA_desc *A,
                            double *W,
                            PLASMA_desc *T,
                            PLASMA_Complex64_t *Q, int LDQ,
                            PLASMA_sequence *sequence, PLASMA_request *request)
{
    plasma_context_t *plasma;
    PLASMA_desc descA;
    PLASMA_desc descT;
    PLASMA_desc descQ;
    PLASMA_Complex64_t *AB;
    double *E;
    int N;
    int NB;
    int LDAB;
    int status;

    plasma = plasma_context_self();
    if (plasma == NULL) {
        plasma_fatal_error("PLASMA_zheevd_Tile_Async", "PLASMA not initialized");
        return PLASMA_ERR_NOT_INITIALIZED;
    }
    if (sequence == NULL) {
        plasma_fatal_error("PLASMA_zheevd_Tile_Async", "NULL sequence");
        return PLASMA_ERR_UNALLOCATED;
    }
    if (request == NULL) {
        plasma_fatal_error("PLASMA_zheevd_Tile_Async", "NULL request");
        return PLASMA_ERR_UNALLOCATED;
    }
    /* Check sequence status */
    if (sequence->status == PLASMA_SUCCESS)
        request->status = PLASMA_SUCCESS;
    else
        return plasma_request_fail(sequence, request, PLASMA_ERR_SEQUENCE_FLUSHED);

    /* Check descriptors for correctness */
    if (plasma_desc_check(A) != PLASMA_SUCCESS) {
        plasma_error("PLASMA_zheevd_Tile_Async", "invalid descriptor");
        return plasma_request_fail(sequence, request, PLASMA_ERR_ILLEGAL_VALUE);
    } else {
        descA = *A;
    }
    if (plasma_desc_check(T) != PLASMA_SUCCESS) {
        plasma_error("PLASMA_zheevd_Tile_Async", "invalid descriptor");
        return plasma_request_fail(sequence, request, PLASMA_ERR_ILLEGAL_VALUE);
    } else {
        descT = *T;
    }
    /* Check input arguments */
    if (jobz != PlasmaNoVec && jobz != PlasmaVec) {
        plasma_error("PLASMA_zheevd_Tile_Async", "illegal value of jobz");
        return plasma_request_fail(sequence, request, PLASMA_ERR_ILLEGAL_VALUE);
    }
    if (uplo != PlasmaLower && uplo != PlasmaUpper) {
        plasma_error("PLASMA_zheevd_Tile_Async", "illegal value of uplo");
        return plasma_request_fail(sequence, request, PLASMA_ERR_ILLEGAL_VALUE);
    }
    if (descA.m != descA.n) {
        plasma_error("PLASMA_zheevd_Tile_Async", "matrix need to be square");
        return plasma_request_fail(sequence, request, PLASMA_ERR_ILLEGAL_VALUE);
    }
    if (descA.nb != descA.mb) {
        plasma_error("PLASMA_zheevd_Tile_Async", "only square tiles supported");
        return plasma_request_fail(sequence, request, PLASMA_ERR_ILLEGAL_VALUE);
    }

    N    = descA.m;
    NB = min(descA.mb,descA.m);
    LDAB = 2*NB+1;

    /* Allocate workspace for band storage of the band matrix A and for the off diagonal after tridiagonalisation */
    AB = (PLASMA_Complex64_t *)plasma_shared_alloc(plasma, LDAB*N, PlasmaComplexDouble);
    memset( AB, 0, LDAB * N * sizeof(PLASMA_Complex64_t) );
    if (AB == NULL) {
        plasma_error("PLASMA_zheevd_Tile_Async", "plasma_shared_alloc(AB) failed");
        plasma_shared_free(plasma, AB);
        return PLASMA_ERR_OUT_OF_RESOURCES;
    }
    E = (double *)plasma_shared_alloc(plasma, N, PlasmaRealDouble);
    if (E == NULL) {
        plasma_error("PLASMA_zheevd_Tile_Async", "plasma_shared_alloc(E) failed");
        plasma_shared_free(plasma, E);
        return PLASMA_ERR_OUT_OF_RESOURCES;
    }
    #if defined(ENABLE_TIMER)
    PLASMA_Double_t timelpk=0.0,timeaplQ2=0.0, timeT=0.0;
    PLASMA_Double_t timeB=0.0,timeblg=0.0,timeaplQ1=0.0,timeconv1=0.0,timeconv2=0.0,timeall=0.0;
    timeall = PLASMA_Wtime();
    timeB   = PLASMA_Wtime();
    #endif
    /* Reduction to tridiagonal form
     * with a two-stage approach.
     */
    /*=======================================
     *  calling Reduction to BAND
     *  then convert matrix to band form
     *=======================================*/
    plasma_dynamic_call_5(plasma_pzhetrd_he2hb,
        PLASMA_enum, uplo,
        PLASMA_desc, descA,
        PLASMA_desc, descT,
        PLASMA_sequence*, sequence,
        PLASMA_request*, request);
    plasma_dynamic_call_6( plasma_pzhbcpy_t2bl,
        PLASMA_enum, uplo,
        PLASMA_desc, descA,
        PLASMA_Complex64_t*, AB,
        int, LDAB,
        PLASMA_sequence*, sequence,
        PLASMA_request*, request);

    plasma_dynamic_sync();

    status = sequence->status;
    if (status != PLASMA_SUCCESS) {
        plasma_error("PLASMA_zheevd","pzhetrd_he2hb+pzcopy");
        return status;
    }
    /*=======================================
     *  END of calling Reduction to BAND
     *=======================================*/
    #if defined(ENABLE_TIMER)
    timeB   = PLASMA_Wtime()-timeB;
    printf("\n  Finish Band red    timing= %lf \n",timeB);
    #endif
    /*=======================================
     *  calling bulge chasing
     *=======================================*/
    PLASMA_Complex64_t *TAU2 = NULL;
    PLASMA_Complex64_t *V2 = NULL;
    PLASMA_Complex64_t *T2 = NULL;
    int Vblksiz, blkcnt, LDT, LDV;
    int WANTZ   = 0;
    int blguplo = PlasmaLower;
    /* int NE      = N; // for later use when a portion of the eigenvectors are requested*/
    if( jobz == PlasmaNoVec )
        WANTZ=0;
    else
        WANTZ=2;
    /* Vblksiz correspond to the blocking used when applying V2 to the matrix Q
     * it is similar to IB in LAPACK ZUNMQR.
     * blkcnt is the number of losange or tile of Vs */
    /* Note that in case PlamaVec requested, the V2 and T2 are stored by the
     * bulgechasing function in a special format:
     * for V2s: it store the V2(LDV,Vblksiz) of each losange in a tile storage meaning
     * that V2_1 is stored then V2_2,..., V2_blkcnt.
     * blkcnt is the number of losange.
     * */
    Vblksiz = min(NB,48);
    LDT     = Vblksiz;
    if( jobz == PlasmaVec ) {
        findVTsiz(N, NB, Vblksiz, &blkcnt, &LDV);
        TAU2   = (PLASMA_Complex64_t *) plasma_shared_alloc(plasma,     blkcnt*Vblksiz, PlasmaComplexDouble);
        V2     = (PLASMA_Complex64_t *) plasma_shared_alloc(plasma, LDV*blkcnt*Vblksiz, PlasmaComplexDouble);
        T2     = (PLASMA_Complex64_t *) plasma_shared_alloc(plasma, LDT*blkcnt*Vblksiz, PlasmaComplexDouble);
        if ( (TAU2 == NULL) || (V2 == NULL) || (T2 == NULL) ) {
            plasma_error("PLASMA_zheevd", "plasma_shared_alloc() failed");
            plasma_shared_free(plasma, TAU2);
            plasma_shared_free(plasma, V2);
            plasma_shared_free(plasma, T2);
            return PLASMA_ERR_OUT_OF_RESOURCES;
        }
        memset(TAU2, 0,     blkcnt*Vblksiz*sizeof(PLASMA_Complex64_t));
        memset(V2,   0, LDV*blkcnt*Vblksiz*sizeof(PLASMA_Complex64_t));
        memset(T2,   0, LDT*blkcnt*Vblksiz*sizeof(PLASMA_Complex64_t));
    }
    else {
        TAU2   = (PLASMA_Complex64_t *) plasma_shared_alloc(plasma, 2*N, PlasmaComplexDouble);
        V2     = (PLASMA_Complex64_t *) plasma_shared_alloc(plasma, 2*N, PlasmaComplexDouble);
        if ( (TAU2 == NULL) || (V2 == NULL) ) {
            plasma_error("PLASMA_zheevd", "plasma_shared_alloc() failed");
            plasma_shared_free(plasma, TAU2);
            plasma_shared_free(plasma, V2);
            return PLASMA_ERR_OUT_OF_RESOURCES;
        }
        memset(TAU2, 0, 2*N*sizeof(PLASMA_Complex64_t));
        memset(V2,   0, 2*N*sizeof(PLASMA_Complex64_t));
    }
    #if defined(ENABLE_TIMER)
    timeblg  = PLASMA_Wtime();
    #endif
    plasma_static_call_13(plasma_pzhetrd_hb2st_v1,
        PLASMA_enum,         blguplo,
        int,                 N,
        int,                 NB,
        int,                 Vblksiz,
        PLASMA_Complex64_t*, AB,
        int,                 LDAB,
        PLASMA_Complex64_t*, V2,
        PLASMA_Complex64_t*, TAU2,
        double*,             W,
        double*,             E,
        int,                 WANTZ,
        PLASMA_sequence*,    sequence,
        PLASMA_request*,     request);
    /* WARNING: If plasma_pzhetrd_hb2st is implemented through a dynamic call, don't
     * forget to synchronize */
    plasma_dynamic_sync();
    /*=======================================
     *  END of calling bulge chasing
     *=======================================*/
    #if defined(ENABLE_TIMER)
    timeblg   = PLASMA_Wtime()-timeblg;
    printf("  Finish Bulge       timing= %lf \n" ,timeblg);
    timelpk   = PLASMA_Wtime();
    #endif
    /*=======================================
     *  calling eigensolver
     *=======================================*/
    /* call eigensolver using lapack routine for our resulting tridiag [D E] */
    plasma_setlapack_multithreads(plasma->world_size);
    /*
    status = LAPACKE_zstevd( LAPACK_COL_MAJOR, lapack_const(jobz), N, W, E, Q, LDQ );
    */
    if(jobz == PlasmaNoVec){
        status = LAPACKE_zstedc( LAPACK_COL_MAJOR, lapack_const(PlasmaNoVec),
                        N, W, E, Q, LDQ );
    } else {
        status = LAPACKE_zstedc( LAPACK_COL_MAJOR, lapack_const(PlasmaIvec),
                        N, W, E, Q, LDQ );
    }

    if(status != 0){
        plasma_error("PLASMA_zstedc","ZSTEDC");
        return status;
    }
    sequence->status = status;
    plasma_setlapack_sequential(plasma);
    /*=======================================
     *  END of calling eigensolver
     *=======================================*/
    #if defined(ENABLE_TIMER)
    timelpk = PLASMA_Wtime()-timelpk;
    printf("  Finish Eigensolver timing= %lf  WANTZ %d threads %d\n" ,timelpk, WANTZ,  plasma->world_size);
    #endif
    if (jobz == PlasmaVec){
        /*=======================================
         *  apply Q2 from the bulge
         *=======================================*/
        /* compute T2 */
        #if defined(ENABLE_TIMER)
        timeT  = PLASMA_Wtime();
        #endif
        plasma_static_call_8(plasma_pzlarft_blgtrd,
            int,                 N,
            int,                 NB,
            int,                 Vblksiz,
            PLASMA_Complex64_t*, V2,
            PLASMA_Complex64_t*, T2,
            PLASMA_Complex64_t*, TAU2,
            PLASMA_sequence*,    sequence,
            PLASMA_request*,     request);
        #if defined(ENABLE_TIMER)
        plasma_dynamic_sync();
        timeT  = PLASMA_Wtime()-timeT;
        printf("  Finish compute T2  timing= %lf \n" ,timeT);
        timeaplQ2  = PLASMA_Wtime();
        #endif
        /* apply Q2 from Left */
        plasma_static_call_14(plasma_pzunmqr_blgtrd,
            PLASMA_enum,         PlasmaLeft,
            PLASMA_enum,         PlasmaNoTrans,
            int,                 N,
            int,                 NB,
            int,                 N,
            int,                 Vblksiz,
            int,                 WANTZ,
            PLASMA_Complex64_t*, V2,
            PLASMA_Complex64_t*, T2,
            PLASMA_Complex64_t*, TAU2,
            PLASMA_Complex64_t*, Q,
            int,                 LDQ,
            PLASMA_sequence*,    sequence,
            PLASMA_request*,     request);
        #if defined(ENABLE_TIMER)
        plasma_dynamic_sync();
        timeaplQ2  = PLASMA_Wtime()-timeaplQ2;
        printf("  Finish compute Q2  timing= %lf \n" ,timeaplQ2);
        #endif
        /*=======================================
         *  apply Q1 from the reduction to band
         *=======================================*/
        /* CASE NB>N, Q1 doesn't need to be applied, only bulge chasing has been done */
        if( NB < N ){
            #if defined(ENABLE_TIMER)
            plasma_dynamic_sync();
            timeconv1   = PLASMA_Wtime();
            #endif
            if ( PLASMA_TRANSLATION == PLASMA_OUTOFPLACE ) {
                plasma_zooplap2tile( descQ, Q, NB, NB, LDQ, N, 0, 0, N, N, sequence, request, plasma_desc_mat_free(&(descQ)) );
            }else {
                plasma_ziplap2tile( descQ, Q, NB, NB, LDQ, N, 0, 0, N, N, sequence, request );
            }
            #if defined(ENABLE_TIMER)
            timeconv1    = PLASMA_Wtime()-timeconv1;
            timeaplQ1   = PLASMA_Wtime();
            #endif
            /* Accumulate the transformations from the first stage */
            if(uplo==PlasmaLower){
                plasma_dynamic_call_7(plasma_pzunmqr,
                    PLASMA_enum, PlasmaLeft,
                    PLASMA_enum, PlasmaNoTrans,
                    PLASMA_desc, plasma_desc_submatrix(descA, descA.mb, 0, descA.m-descA.mb, descA.n-descA.nb),
                    PLASMA_desc, plasma_desc_submatrix(descQ, descQ.mb, 0, descQ.m-descQ.mb, descQ.n),
                    PLASMA_desc, plasma_desc_submatrix(descT, descT.mb, 0, descT.m-descT.mb, descT.n-descT.nb),
                    PLASMA_sequence*, sequence,
                    PLASMA_request*, request);

            }
            else {
                plasma_dynamic_call_7(plasma_pzunmlq,
                    PLASMA_enum, PlasmaLeft,
                    PLASMA_enum, PlasmaConjTrans,
                    PLASMA_desc, plasma_desc_submatrix(descA, 0, descA.nb, descA.m-descA.mb, descA.n-descA.nb),
                    PLASMA_desc, plasma_desc_submatrix(descQ, descQ.mb, 0, descQ.m-descQ.mb, descQ.n         ),
                    PLASMA_desc, plasma_desc_submatrix(descT, 0, descT.nb, descT.m-descT.mb, descT.n-descT.nb),
                    PLASMA_sequence*, sequence,
                    PLASMA_request*, request);
            }
            #if defined(ENABLE_TIMER)
            plasma_dynamic_sync();
            timeaplQ1   = PLASMA_Wtime()-timeaplQ1;
            printf("  Finish compute Q1  timing= %lf \n" ,timeaplQ1);
            timeconv2    = PLASMA_Wtime();
            #endif
            if ( PLASMA_TRANSLATION == PLASMA_OUTOFPLACE ) {
                plasma_zooptile2lap( descQ, Q, NB, NB, LDQ, N, sequence, request );
                plasma_dynamic_sync();
                plasma_desc_mat_free(&descQ);
            } else {
                plasma_ziptile2lap( descQ, Q, NB, NB, LDQ, N, sequence, request );
                plasma_dynamic_sync();
            }
            #if defined(ENABLE_TIMER)
            plasma_dynamic_sync();
            timeconv2    = PLASMA_Wtime()-timeconv2;
            printf("  Finish convert     timing= %lf \n" ,timeconv1+timeconv2);
            #endif
        } /* END of ( NB < N ) */
    }
    #if defined(ENABLE_TIMER)
    timeall = PLASMA_Wtime()-timeall;
    printf("  Finish full eigenproblem threads %d  N %d  timeall= %lf \n", plasma->world_size, N, timeall);
    #endif

    if( jobz == PlasmaVec )
        plasma_shared_free(plasma, T2);
    plasma_shared_free(plasma, V2);
    plasma_shared_free(plasma, TAU2);
    plasma_shared_free(plasma, E);
    plasma_shared_free(plasma, AB);
    return PLASMA_SUCCESS;
}

ext: c


text:
/* See compilation notes in simple_cplx.c*/
#include "cplx.h"                 //gsl_cplx_from_c99; see below.
#include <gsl/gsl_blas.h>         //gsl_blas_ddot
#include <gsl/gsl_complex_math.h> //gsl_complex_mul(_real)

gsl_vector_complex *cvec_dot_gslcplx(gsl_vector_complex *v, gsl_complex x){
    gsl_vector_complex *out = gsl_vector_complex_alloc(v->size);
    for (int i=0; i< v->size; i++)
        gsl_vector_complex_set(out, i,
                            gsl_complex_mul(x, gsl_vector_complex_get(v, i)));
    return out;
}

gsl_vector_complex *vec_dot_gslcplx(gsl_vector *v, gsl_complex x){
    gsl_vector_complex *out = gsl_vector_complex_alloc(v->size);
    for (int i=0; i< v->size; i++)
        gsl_vector_complex_set(out, i,
                            gsl_complex_mul_real(x, gsl_vector_get(v, i)));
    return out;
}

gsl_vector_complex *cvec_dot_c(gsl_vector_complex *v, complex double x){
    return cvec_dot_gslcplx(v, gsl_cplx_from_c99(x));
}

gsl_vector_complex *vec_dot_c(gsl_vector *v, complex double x){
    return vec_dot_gslcplx(v, gsl_cplx_from_c99(x));
}

complex double ddot (complex double x, complex double y){return x*y;}

void gsl_vector_complex_print(gsl_vector_complex *v){
    for (int i=0; i< v->size; i++) {
        gsl_complex x = gsl_vector_complex_get(v, i);
        printf("%4g+%4gi%c", GSL_REAL(x), GSL_IMAG(x), i < v->size-1 ? '\t' : '\n');
    }
}

max_stars_repo_name: catalingheorghe/learning-activities

max_issues_repo_name: catalingheorghe/learning-activities

max_forks_repo_name: catalingheorghe/learning-activities


text:
/**
 * The interface to the binning and weight routines
 */
#ifndef _SPCE_BINNING_H
#define _SPCE_BINNING_H 1

#include <math.h>
#include <gsl/gsl_vector.h>

#include "spce_output.h"
#include "spc_trace_functions.h"
#include "aXe_grism.h"
#include "spc_spc.h"

/**
  defines a pixel for weighting purposes.  Here, p0, p1, p2, p3 are
  the coordinates of the points where the square first contribute,
  where its contribution becomes constant, where it slopes down again,
  and where it doesn't contribute any more.  For example, a pixel at (0,0)
  viewed at from a side has p0=p1=0, p2=p3=1, whereas the same pixel
  viewed from a corner has p0=-1, p1=p2=0, p3=1;  slope is the
  ascent on the non-constant parts, fmax the maximum contribution.
  angle is the angle viewing.  Its tangent tana or its cotangent
  cota may be undefined for certain viewing angles.  The function
  weight_function is selected such that this is unimportatnt.
  x0, y0, and size are the coordinates of the lower left corner and
  the size of the square.
*/
typedef struct w_pixel_s
{
     double p0, p1, p2, p3;
     double angle;
     double tana, cota;
     double x0, y0, size;
     double fmax, slope;
     double (*weight_function) (const double x1, const double y1,
				const double x2, const double y2,
				const struct w_pixel_s * const pix);
}
w_pixel;

extern spectrum *
bin_naive (const ap_pixel * const ap_p, const double ob_width,
	   const double ob_orient, const int quant_cont);

extern spectrum *
bin_weighted (const ap_pixel * const ap_p, const double ob_orient,
	      const trace_func * const trace, const int n_sub,
	      const int flags);

extern spectrum *
bin_optimal (const ap_pixel * const ap_p, const beam curbeam,
	     const int quant_cont, const gsl_matrix *weights,
	     const drzstamp_dim dimension, gsl_matrix *coverage);

#endif

max_stars_repo_name: sosey/pyaxe

max_issues_repo_name: sosey/pyaxe

max_forks_repo_name: sosey/pyaxe


text:
#include "simple_blas.h"

#ifdef ACCELERATE
#include <Accelerate/Accelerate.h>
#else
#include <cblas.h>
#include "lapack.h"
#endif

int to_cblas_transpose(const morpheus_transpose_e trans) {
  switch (trans) {
  case morpheus_no_trans:
    return CblasNoTrans;
  case morpheus_trans:
    return CblasTrans;
  }
  return 0;
}

int to_cblas_layout(const morpheus_layout_e layout) {
  switch (layout) {
  case morpheus_row_major:
    return CblasRowMajor;
  case morpheus_col_major:
    return CblasColMajor;
  }
  return 0;
}

void morpheus_dscal(const int n, const double alpha, double *x) {
  const int inc = 1;
  cblas_dscal(n, alpha, x, inc);
}

void morpheus_dcopy(const int n, const double *x, double *y) {
  const int inc = 1;
  cblas_dcopy(n, x, inc, y, inc);
}

void morpheus_dcopy_scalar(const int n, const double alpha, double *x) {
  cblas_dcopy(n, &alpha, 0, x, 1);
}

void morpheus_dgemv(const morpheus_layout_e layout,
                    const morpheus_transpose_e trans,
                    const int nrows, const int ncols,
                    const double alpha, const double *A, const double *x,
                    const double beta, double *y) {
  const int inc = 1;
  const int lda = layout == morpheus_row_major ? ncols : nrows;
  cblas_dgemv(to_cblas_layout(layout), to_cblas_transpose(trans), nrows, ncols,
              alpha, A, lda, x, inc, beta, y, inc);
}

void morpheus_daxpy(const int n, const double alpha,
                    const double *x, double *y) {
  const int inc = 1;
  cblas_daxpy(n, alpha, x, inc, y, inc);
}

double morpheus_dnrm2(const int n, const double *x) {
  const int inc = 1;
  return cblas_dnrm2(n, x, inc);
}

double morpheus_ddot(const int n, const double *x, const double *y) {
  const int inc = 1;
  return cblas_ddot(n, x, inc, y, inc);
}

int morpheus_inverse(int n, double *x, int *pivot, double *workspace) {
  workspace = workspace+n;
  int lwork = n*n;

  int rc;
  dgetrf_(&n, &n, x, &n, pivot, &rc);
  if (rc != 0) {
    return rc;
  }
  dgetri_(&n, x, &n, pivot, workspace, &lwork, &rc);
  return rc;
}

void morpheus_dger(const morpheus_layout_e layout,
                   const int nrows,
                   const int ncols,
                   const double alpha,
                   const double *x,
                   const double *y,
                   double *a) {
  const int inc = 1;
  const int lda = layout == morpheus_row_major ? ncols : nrows;
  cblas_dger(to_cblas_layout(layout), nrows, ncols,
                             alpha, x, inc, y, inc, a, lda);
}

void morpheus_dgemm(const morpheus_layout_e layout,
                    const morpheus_transpose_e trans_a,
                    const morpheus_transpose_e trans_b,
                    const int m,
                    const int n,
                    const int k,
                    const double alpha,
                    const double *a,
                    const double *b,
                    const double beta,
                    double *c) {
  int lda, ldb, ldc;
  if (layout == morpheus_row_major) {
    lda = trans_a == morpheus_no_trans ? k : m;
    ldb = trans_b == morpheus_no_trans ? n : k;
    ldc = n;
  } else {
    lda = trans_a == morpheus_no_trans ? m : k;
    ldb = trans_b == morpheus_no_trans ? k : n;
    ldc = m;
  }
  cblas_dgemm(to_cblas_layout(layout),
              to_cblas_transpose(trans_a), to_cblas_transpose(trans_b),
              m, n, k, alpha, a, lda, b, ldb, beta, c, ldc);
}

void morpheus_identity(const int n, double *a) {
  const int m = n*n;
  const double zero = 0;
  const double one = 1;
  cblas_dcopy(m, &zero, 0, a, 1);
  cblas_dcopy(n, &one, 0, a, n+1);
}

max_stars_repo_name: Alexander-Ignatyev/morpheus

max_issues_repo_name: aligusnet/morpheus

max_forks_repo_name: aligusnet/morpheus


text:
#if !defined(fvSupport_h)
#define fvSupport_h

#include <petsc.h>

/**
 * reproduces the petsc call with grad fixes for multiple fields
 * @param dm
 * @param fvm
 * @param locX
 * @param grad
 * @return
 */
PETSC_EXTERN PetscErrorCode DMPlexGetDataFVM_MulfiField(DM dm, PetscFV fv, Vec *cellgeom, Vec *facegeom, DM *gradDM);

/**
 * Check to make sure local ghost boundary have valid gradient values
 */
PETSC_EXTERN PetscErrorCode ABLATE_FillGradientBoundary(DM dm, PetscFV auxFvm, Vec localXVec, Vec gradLocalVec);

#endif

max_stars_repo_name: mtmcgurn-buffalo/ablate

max_issues_repo_name: mtmcgurn-buffalo/ablate

max_forks_repo_name: mtmcgurn-buffalo/ablate


text:
/* interpolation/bilinear.c
 *
 * Copyright 2012 David Zaslavsky
 *
 * 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.
 */

#include <config.h>
#include <gsl/gsl_errno.h>
#include <gsl/gsl_interp.h>
#include <gsl/gsl_interp2d.h>

#define IDX2D(i, j, xsize, ysize) ((j) * (xsize) + (i))

static int
bilinear_init(void * state, const double xa[], const double ya[],
              const double za[], size_t xsize, size_t ysize)
{
  return GSL_SUCCESS;
}

static int
bilinear_eval(const void * state, const double xarr[], const double yarr[],
              const double zarr[], size_t xsize, size_t ysize,
              double x, double y, gsl_interp_accel * xa,
              gsl_interp_accel * ya, double * z)
{
  double xmin, xmax, ymin, ymax, zminmin, zminmax, zmaxmin, zmaxmax;
  double dx, dy;
  double t, u;
  size_t xi, yi;

  if (xa != NULL)
    xi = gsl_interp_accel_find(xa, xarr, xsize, x);
  else
    xi = gsl_interp_bsearch(xarr, x, 0, xsize - 1);

  if (ya != NULL)
    yi = gsl_interp_accel_find(ya, yarr, ysize, y);
  else
    yi = gsl_interp_bsearch(yarr, y, 0, ysize - 1);

  xmin = xarr[xi];
  xmax = xarr[xi + 1];
  ymin = yarr[yi];
  ymax = yarr[yi + 1];
  zminmin = zarr[IDX2D(xi, yi, xsize, ysize)];
  zminmax = zarr[IDX2D(xi, yi + 1, xsize, ysize)];
  zmaxmin = zarr[IDX2D(xi + 1, yi, xsize, ysize)];
  zmaxmax = zarr[IDX2D(xi + 1, yi + 1, xsize, ysize)];
  dx = xmax - xmin;
  dy = ymax - ymin;
  t = (x - xmin)/dx;
  u = (y - ymin)/dy;
  *z = (1.-t)*(1.-u)*zminmin + t*(1.-u)*zmaxmin + (1.-t)*u*zminmax + t*u*zmaxmax;

  return GSL_SUCCESS;
}

static int
bilinear_deriv_x(const void * state, const double xarr[],
                 const double yarr[], const double zarr[],
                 size_t xsize, size_t ysize, double x, double y,
                 gsl_interp_accel * xa, gsl_interp_accel * ya, double * z_p)
{
  double xmin, xmax, ymin, ymax, zminmin, zminmax, zmaxmin, zmaxmax;
  double dx, dy;
  double dt, u;
  size_t xi, yi;

  if (xa != NULL)
    xi = gsl_interp_accel_find(xa, xarr, xsize, x);
  else
    xi = gsl_interp_bsearch(xarr, x, 0, xsize - 1);

  if (ya != NULL)
    yi = gsl_interp_accel_find(ya, yarr, ysize, y);
  else
    yi = gsl_interp_bsearch(yarr, y, 0, ysize - 1);

  xmin = xarr[xi];
  xmax = xarr[xi + 1];
  ymin = yarr[yi];
  ymax = yarr[yi + 1];
  zminmin = zarr[IDX2D(xi, yi, xsize, ysize)];
  zminmax = zarr[IDX2D(xi, yi + 1, xsize, ysize)];
  zmaxmin = zarr[IDX2D(xi + 1, yi, xsize, ysize)];
  zmaxmax = zarr[IDX2D(xi + 1, yi + 1, xsize, ysize)];
  dx = xmax - xmin;
  dy = ymax - ymin;
  dt = 1./dx; /* partial t / partial x */
  u = (y - ymin)/dy;
  *z_p = dt*(-(1.-u)*zminmin + (1.-u)*zmaxmin - u*zminmax + u*zmaxmax);

  return GSL_SUCCESS;
}

static int
bilinear_deriv_y(const void * state, const double xarr[],
                 const double yarr[], const double zarr[],
                 size_t xsize, size_t ysize, double x, double y,
                 gsl_interp_accel * xa, gsl_interp_accel * ya, double * z_p)
{
  double xmin, xmax, ymin, ymax, zminmin, zminmax, zmaxmin, zmaxmax;
  double dx, dy;
  double t, du;
  size_t xi, yi;

  if (xa != NULL)
    xi = gsl_interp_accel_find(xa, xarr, xsize, x);
  else
    xi = gsl_interp_bsearch(xarr, x, 0, xsize - 1);

  if (ya != NULL)
    yi = gsl_interp_accel_find(ya, yarr, ysize, y);
  else
    yi = gsl_interp_bsearch(yarr, y, 0, ysize - 1);

  xmin = xarr[xi];
  xmax = xarr[xi + 1];
  ymin = yarr[yi];
  ymax = yarr[yi + 1];
  zminmin = zarr[IDX2D(xi, yi, xsize, ysize)];
  zminmax = zarr[IDX2D(xi, yi + 1, xsize, ysize)];
  zmaxmin = zarr[IDX2D(xi + 1, yi, xsize, ysize)];
  zmaxmax = zarr[IDX2D(xi + 1, yi + 1, xsize, ysize)];
  dx = xmax - xmin;
  dy = ymax - ymin;
  t = (x - xmin)/dx;
  du = 1./dy; /* partial u / partial y */
  *z_p = du*(-(1.-t)*zminmin - t*zmaxmin + (1.-t)*zminmax + t*zmaxmax);

  return GSL_SUCCESS;
}

static int
bilinear_deriv2(const void * state, const double xarr[],
                const double yarr[], const double zarr[],
                size_t xsize, size_t ysize, double x, double y,
                gsl_interp_accel * xa, gsl_interp_accel * ya, double * z_pp)
{
  *z_pp = 0.0;
  return GSL_SUCCESS;
}

static int
bilinear_derivxy(const void * state, const double xarr[],
                 const double yarr[], const double zarr[],
                 size_t xsize, size_t ysize, double x, double y,
                 gsl_interp_accel * xa, gsl_interp_accel * ya, double * z_pp)
{
  double xmin, xmax, ymin, ymax, zminmin, zminmax, zmaxmin, zmaxmax;
  double dx, dy;
  double dt, du;
  size_t xi, yi;

  if (xa != NULL)
    xi = gsl_interp_accel_find(xa, xarr, xsize, x);
  else
    xi = gsl_interp_bsearch(xarr, x, 0, xsize - 1);

  if (ya != NULL)
    yi = gsl_interp_accel_find(ya, yarr, ysize, y);
  else
    yi = gsl_interp_bsearch(yarr, y, 0, ysize - 1);

  xmin = xarr[xi];
  xmax = xarr[xi + 1];
  ymin = yarr[yi];
  ymax = yarr[yi + 1];
  zminmin = zarr[IDX2D(xi, yi, xsize, ysize)];
  zminmax = zarr[IDX2D(xi, yi + 1, xsize, ysize)];
  zmaxmin = zarr[IDX2D(xi + 1, yi, xsize, ysize)];
  zmaxmax = zarr[IDX2D(xi + 1, yi + 1, xsize, ysize)];
  dx = xmax - xmin;
  dy = ymax - ymin;
  dt = 1./dx; /* partial t / partial x */
  du = 1./dy; /* partial u / partial y */
  *z_pp = dt*du*(zminmin-zmaxmin-zminmax+zmaxmax);

  return GSL_SUCCESS;
}

static const gsl_interp2d_type bilinear_type = {
  "bilinear",
  2,
  NULL,
  &bilinear_init,
  &bilinear_eval,
  &bilinear_deriv_x,
  &bilinear_deriv_y,
  &bilinear_deriv2,
  &bilinear_derivxy,
  &bilinear_deriv2,
  NULL
};

const gsl_interp2d_type * gsl_interp2d_bilinear = &bilinear_type;

#undef IDX2D

max_stars_repo_name: ielomariala/Hex-Game

max_issues_repo_name: Brian-ning/HMNE

max_forks_repo_name: Brian-ning/HMNE


text:
#ifndef PETSC4PY_COMPAT_H
#define PETSC4PY_COMPAT_H

#include <petsc.h>
#include "compat/mpi.h"
#include "compat/hdf5.h"
#include "compat/mumps.h"
#include "compat/hypre.h"
#include "compat/tao.h"

#endif/*PETSC4PY_COMPAT_H*/