text:
#pragma once
#include "errors.h"
#include "usageinfoformat.h"
#include "detail/config.h"
#include "detail/configmacro.h"
#include "detail/configaccess.h"
#include "detail/flag.h"
#include <gsl/gsl>
#include <iostream>
#include <ostream>
#include <optional>
#include <map>
#include <utility>

namespace cmdlime{

enum class ErrorOutputMode{
    STDOUT,
    STDERR
};

template<typename TConfig>
class ConfigReader{
    struct CommandHelpFlag{
        bool value = false;
        std::string usageInfo;
    };

public:
    ConfigReader(TConfig& cfg,
                 std::string programName,
                 const UsageInfoFormat& usageInfoFormat = {},
                 ErrorOutputMode errorOutputMode = ErrorOutputMode::STDERR)
        : cfg_(cfg)
        , programName_(std::move(programName))
        , usageInfoFormat_(usageInfoFormat)
        , errorOutput_(errorOutputMode == ErrorOutputMode::STDERR ? std::cerr : std::cout)
    {}

    int exitCode() const
    {
        return exitCode_;
    }

    bool read(const std::vector<std::string>& cmdLine)
    {
        addExitFlags();
        if(!processCommandLine(cmdLine))
            return exitOnError(-1);

        if (processFlagsAndExit())
            return exitOnFlag();

        return success();
    }

    bool readCommandLine(int argc, char** argv)
    {
        auto cmdLine = std::vector<std::string>(argv + 1, argv + argc);
        return read(cmdLine);
    }

private:
    void addExitFlags()
    {
        using NameProvider = typename detail::Format<detail::ConfigAccess<TConfig>::format()>::nameProvider;
        auto helpFlag = std::make_unique<detail::Flag>(NameProvider::name("help"),
                                                       std::string{},
                                                       [this]()->bool&{return help_;},
                                                       detail::Flag::Type::Exit);
        helpFlag->info().addDescription("show usage info and exit");
        detail::ConfigAccess<TConfig>{cfg_}.addFlag(std::move(helpFlag));

        if (!cfg_.versionInfo().empty()){
            auto versionFlag = std::make_unique<detail::Flag>(NameProvider::name("version"),
                                                           std::string{},
                                                           [this]()->bool&{return version_;},
                                                           detail::Flag::Type::Exit);
            versionFlag->info().addDescription("show version info and exit");
            detail::ConfigAccess<TConfig>{cfg_}.addFlag(std::move(versionFlag));
        }

        detail::ConfigAccess<TConfig>(cfg_).addHelpFlagToCommands(programName_);
    }

    bool processCommandLine(const std::vector<std::string>& cmdLine)
    {
        try{
            cfg_.read(cmdLine);
        }
        catch(const CommandError& e){
            errorOutput_ << "Command '" + e.commandName() + "' error: " << e.what() << "\n";
            std::cout << e.commandUsageInfo() << std::endl;
            return false;
        }
        catch(const Error& e){
            errorOutput_ << e.what() << "\n";
            std::cout << cfg_.usageInfo(programName_) << std::endl;
            return false;
        }
        return true;
    }

    bool processFlagsAndExit()
    {
        if (help_){
            std::cout << cfg_.usageInfoDetailed(programName_, usageInfoFormat_) << std::endl;
            return true;
        }
        if (version_){
            std::cout << cfg_.versionInfo() << std::endl;
            return true;
        }

        for (auto command :  detail::ConfigAccess<TConfig>{cfg_}.commandList())
            if (checkCommandHelpFlag(command))
                return true;

        return false;
    }

    bool checkCommandHelpFlag(gsl::not_null<detail::ICommand*> command)
    {
        if (command->isHelpFlagSet()){
            std::cout << command->usageInfoDetailed() << std::endl;
            return true;
        }

        for (auto childCommand : command->commandList())
            if (checkCommandHelpFlag(childCommand))
                return true;

        return false;
    }

    bool exitOnFlag()
    {
        exitCode_ = 0;
        return false;
    }

    bool exitOnError(int errorCode)
    {
        exitCode_ = errorCode;
        return false;
    }

    bool success()
    {
        exitCode_ = 0;
        return true;
    }

private:
    TConfig& cfg_;
    std::string programName_;
    UsageInfoFormat usageInfoFormat_;
    std::map<int, CommandHelpFlag> commandHelpFlags_;
    std::ostream& errorOutput_;
    int exitCode_ = 0;
    bool help_ = false;
    bool version_ = false;
}
;
}

ext: h

max_stars_repo_name: GerHobbelt/hypertextcpp

max_issues_repo_name: GerHobbelt/hypertextcpp

max_forks_repo_name: GerHobbelt/hypertextcpp


text:
#include <math.h>
#include <stdio.h>
#include <stddef.h>
#include <stdlib.h>
#include <string.h>
#include <ctype.h>
#include <gsl/gsl_math.h>
#include <gsl/gsl_blas.h>
#include <gsl/gsl_vector.h>
#include <gsl/gsl_matrix.h>
#include <gsl/gsl_multimin.h>
#include <gsl/gsl_sf.h>
#include <gsl/gsl_randist.h>
#include <gsl/gsl_cdf.h>
#include <gsl/gsl_linalg.h>
#include <getopt.h>
#include <string.h>

/******************************************************************************/
#define H 0.5
#define maxnd 2000

#define neutral 0
#define selected 1
#define max_n_int_evaluated 1000
#define max_n_s_evaluated 2000

#define undefined -99999999999.99
#define undefined_int 999999
#define max_str_len 100
#define n_decimal_places 6
#define max_s_ranges 100

#define s_evaluated_vec_file_const "s_evaluated.dat"
#define t2_evaluated_vec_file_const "t2_evaluated.dat"
#define n2_evaluated_vec_file_const "n2_evaluated.dat"
#define s_range_file_const "s_ranges_evaluated.dat"
#define phase_1_dir_const "n1_n2_t2_s_vec_dir"
#define phase_2_dir_const "n2_t2_s_vec_dir"
#define const_pop_file "control_generate_allele_distrib_n1.out.bin"
/******************************************************************************/
/*Function List*/
double compute_weighted_average(double *weighted_vec, double *phase1_mut_vec,
				double *phase2_mut_vec, int n1d, int n2d);
float assign_float_from_read_buff(char *buffer, int start);
int get_average_egf_vecs_1_and_2(double s,int n1, int n2, double t2_real,
				 int t2_lower, int t2_upper, double *egf_vec, char *buffer_p1_n2_t2_lower,
				 char *buffer_p2_n2_t2_lower, char *buffer_p1_n2_t2_upper,
				 char *buffer_p2_n2_t2_upper);
				 get_upper_lower_int(double parm, int *lower, int *upper, int n_int_evaluated, int *int_evaluated_vec);
get_upper_lower_double(double parm, double*lower, double *upper,int n_int_evaluated, double *double_evaluated_vec);
get_s_ranges();
get_int_evaluated_vec(int *int_evaluated_vec, int *n_int_evaluated,
		      int *int_lower, int *int_step, char *int_evaluated_vec_file);
scale_vector(double *vec, int n2, double total_density_s0);
int get_binary_egf_vec(char *buffer, int n, int t2, double s, double *egf_vec);

fold_vector_double(double *vec, int counter);
fold_vector_int(int *vec, int counter);

int nearest_n2_ind(int n2_real);
read_phase1_phase2_file_into_buffer(int n1, int phase, int n2, int t2, char *buffer,
				    int file_size_bytes);
read_const_pop_file_into_buffer(int n1, char *buffer, int file_size_bytes);
compute_weighted_average_egf_vec(double t2_real, int t2_lower, int t2_upper,
				 double *egf_vec_lower, double *egf_vec_upper,
				 double *egf_vec, int n2d);
double compute_weighted_average(double *weighted_vec, double *phase1_mut_vec,
				double *phase2_mut_vec, int n1d, int n2d);
int get_const_pop_egf_vec(double s, int n1, double *egf_vec, char *buffer_const_pop,
			  int assign_tds0_mode);
assign_int_from_read_buff(int *res, char *buffer, int start);
float assign_float_from_read_buff(char *buffer, int start);
int compute_s_tab_num(double s);
int compute_file_size_bytes(int n);
set_up_file_name(int n1, char *froot, char *file_name);
get_data_path(char *data_path);
get_s_evaluated_vec(double *s_evaluated_vec, int *n_s_evaluated,
   int *n_s_evaluated_file, char *s_evaluated_vec_file);
get_lower_upper(int ind, double *s_evaluated_vec, int n_s_evaluated,
   double *lower, double *upper);
dumpvector(double *v, int min, int max, char *s);
dumpmatrix(double **m, int s1,
        int rows, int s2, int cols, char *s);
double selfn(double q, double s, double h);
setuprow(double **a, double p, int row, int n);
setupmatrix(double **a, int n1, int n2, double s, double h);
matrixinvert(int N, double **a,double *inverter);
tmiterate(double **a,double *mut1,
        int t, int n1, int n2,int decay,double *sumf);
tmiterate2(double **a,double *mut1,
        int t, int n1, int n2,int decay,double **a1);
eqf_using_matrix_inversion(int n1,double s,double **a,double *egf_out);
vector_average(int n1,int n2,double *fv1,double *fv2, double *fv_out);
vector_s_average(int n,double P1,double *fv1,double *fv2, double *fv_out);

output_sfs_to_file_thanasis_format(int n, int *sfs1,int *sfs2,char *sfs_filename);
output_sfs_to_file_peter_format(int n,int sample1, int sample2, int *sfs1,int *sfs2,char *sfs_filename);
get_sfs(int alleles,int *par,char *sfs_filename);

get_sfs_peter1(int *total_neu,int *total_sel,int *nalleles,int *sfs_sel,int *sfs_neu,char *sfs_filename);
get_sfs_peter2(int *total_neu,int *total_sel,int *nalleles,int *sfs_sel,int *sfs_neu,char *sfs_filename,int *sel_sites,int *sel_diff,int *neu_sites,int *neu_diff);

int binarysearchcumulative_from_zero(double *array, int size);
double randomly_select_allele_frequency(double *egf_vec, double *temp1, int nd);
set_up_cumulative_allele_freq_vec(int nd, double *egf_vec, double *cum_vec,
				  char *str);
generate_simulated_data(double *egf_vec, int nsites, int sample_size, int nd,
			int *sfs);

egf_scaling_s(int n1,double f0,double *v_s_in, double *v_s_out);
egf_scaling_f0(int n1,double f0,double *fv);
egf_scaling(int N,double f0,double *FV0,double *FVS);
binomial_sampling(int n1,int alleles,int sample, double *invy,int *discrete);

int set_up_density_vec_equal_effects(double s, int n_s_evaluated, double *gamma_density_vec);
int set_up_density_vec_step_effects(double alpha, double beta, int n_s_evaluated, double *gamma_density_vec);
double compute_beta_densities(double alpha, double *s_evaluated_vec,
			      int n_s_evaluated, double *gamma_density_vec, double beta);
double compute_gamma_densities(double alpha, double *s_evaluated_vec,
				int n_s_evaluated, double *gamma_density_vec, double beta);

double calculate_ne(int n1,int n2,int t);
/******************************************************************************/
/*Structure List*/
struct rangestruct
{
  double s_lower, s_upper, s_step;
  int s_lower_index;
};

struct rangestruct s_ranges[max_s_ranges+1];
/******************************************************************************/
/*Global variable List*/
double allele_freq_sampled[maxnd+1];
double total_density_s0;
char s_evaluated_vec_file[max_str_len], t2_evaluated_vec_file[max_str_len],
  n2_evaluated_vec_file[max_str_len], s_range_file[max_str_len],
  phase_1_dir[max_str_len], phase_2_dir[max_str_len],
  const_pop_dir[max_str_len];
char data_path[max_str_len], file_label_str[max_str_len];
int nranges, trace_level, n_sfs;
int n_s_evaluated, n_s_evaluated_file, gss_n2;
int t2_evaluated_vec[max_n_int_evaluated], n_t2_evaluated, t2_lower;
int t2_step, verbose_mode, neutrals_only_mode;
int n2_evaluated_vec[max_n_int_evaluated], n_n2_evaluated, n2_lower;
double s_evaluated_vec[max_n_s_evaluated];
/******************************************************************************/


text:
// MIT License
//
// Copyright (c) 2020 SunnyCase
//
// Permission is hereby granted, free of charge, to any person obtaining a copy
// of this software and associated documentation files (the "Software"), to deal
// in the Software without restriction, including without limitation the rights
// to use, copy, modify, merge, publish, distribute, sublicense, and/or sell
// copies of the Software, and to permit persons to whom the Software is
// furnished to do so, subject to the following conditions:
//
// The above copyright notice and this permission notice shall be included in all
// copies or substantial portions of the Software.
//
// THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND, EXPRESS OR
// IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES OF MERCHANTABILITY,
// FITNESS FOR A PARTICULAR PURPOSE AND NONINFRINGEMENT. IN NO EVENT SHALL THE
// AUTHORS OR COPYRIGHT HOLDERS BE LIABLE FOR ANY CLAIM, DAMAGES OR OTHER
// LIABILITY, WHETHER IN AN ACTION OF CONTRACT, TORT OR OTHERWISE, ARISING FROM,
// OUT OF OR IN CONNECTION WITH THE SOFTWARE OR THE USE OR OTHER DEALINGS IN THE
// SOFTWARE.
#pragma once
#include <chino/ddk/kernel.h>
#include <cstdint>
#include <gsl/gsl-lite.hpp>
#include <xmmintrin.h>

namespace chino::arch
{
// No asynchronous interrupt is supported by this target
// Only save callee saved registers in context
struct win32_thread_context
{
    uintptr_t rbx;
    uintptr_t rbp;
    uintptr_t rdi;
    uintptr_t rsi;
    uintptr_t rsp;
    uintptr_t r12;
    uintptr_t r13;
    uintptr_t r14;
    uintptr_t r15;
    uintptr_t reserved0;

    __m128 xmm6;
    __m128 xmm7;
    __m128 xmm8;
    __m128 xmm9;
    __m128 xmm10;
    __m128 xmm11;
    __m128 xmm12;
    __m128 xmm13;
    __m128 xmm14;
    __m128 xmm15;

    uintptr_t stack_low;
    uintptr_t stack_high;
};

using thread_context_t = win32_thread_context;

struct win32_arch
{
    static constexpr size_t ALLOCATE_ALIGNMENT = 16;

    static uint32_t current_processor() noexcept { return 0; }
    static void yield_processor() noexcept;

    static uintptr_t disable_irq() noexcept;
    static void restore_irq(uintptr_t state) noexcept;

    static void init_thread_context(thread_context_t &context, gsl::span<uintptr_t> stack, kernel::thread_thunk_t start, void *arg0, void *arg1) noexcept;
    [[noreturn]] static void start_schedule(thread_context_t &context) noexcept;
    static void yield(thread_context_t &old_context, thread_context_t &new_context) noexcept;

    static void init_stack_check() noexcept;
};

using arch_t = win32_arch;
}

max_stars_repo_name: chino-os/chino-os

max_issues_repo_name: dotnetGame/chino-os

max_forks_repo_name: dotnetGame/chino-os


text:
/* ode-initval/gear1.c
 *
 * Copyright (C) 1996, 1997, 1998, 1999, 2000 Gerard Jungman
 *
 * 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.
 */

/* Gear 1. This is the implicit Euler a.k.a backward Euler method. */

/* Author:  G. Jungman
 */

/* Error estimation by step doubling, see eg. Ascher, U.M., Petzold,
   L.R., Computer methods for ordinary differential and
   differential-algebraic equations, SIAM, Philadelphia, 1998.
   The method is also described in eg. this reference.
*/

#include <config.h>
#include <stdlib.h>
#include <string.h>
#include <gsl/gsl_math.h>
#include <gsl/gsl_errno.h>
#include <gsl/gsl_odeiv.h>

#include "odeiv_util.h"

typedef struct
{
  double *k;
  double *y0;
  double *y0_orig;
  double *y_onestep;
}
gear1_state_t;

static void *
gear1_alloc (size_t dim)
{
  gear1_state_t *state = (gear1_state_t *) malloc (sizeof (gear1_state_t));

  if (state == 0)
    {
      GSL_ERROR_NULL ("failed to allocate space for gear1_state", GSL_ENOMEM);
    }

  state->k = (double *) malloc (dim * sizeof (double));

  if (state->k == 0)
    {
      free (state);
      GSL_ERROR_NULL ("failed to allocate space for k", GSL_ENOMEM);
    }

  state->y0 = (double *) malloc (dim * sizeof (double));

  if (state->y0 == 0)
    {
      free (state->k);
      free (state);
      GSL_ERROR_NULL ("failed to allocate space for y0", GSL_ENOMEM);
    }

  state->y0_orig = (double *) malloc (dim * sizeof (double));

  if (state->y0_orig == 0)
    {
      free (state->y0);
      free (state->k);
      free (state);
      GSL_ERROR_NULL ("failed to allocate space for y0_orig", GSL_ENOMEM);
    }

  state->y_onestep = (double *) malloc (dim * sizeof (double));

  if (state->y_onestep == 0)
    {
      free (state->y0_orig);
      free (state->y0);
      free (state->k);
      free (state);
      GSL_ERROR_NULL ("failed to allocate space for y_onestep", GSL_ENOMEM);
    }

  return state;
}

static int
gear1_step (double *y, gear1_state_t *state,
	    const double h, const double t,
	    const size_t dim, const gsl_odeiv_system *sys)
{
  /* Makes an implicit Euler advance with step size h.
     y0 is the initial values of variables y.

     The implicit matrix equations to solve are:

     k = y0 + h * f(t + h, k)

     y = y0 + h * f(t + h, k)
  */

  const int iter_steps = 3;
  int nu;
  size_t i;
  double *y0 = state->y0;
  double *k = state->k;

  /* Iterative solution of k = y0 + h * f(t + h, k)

     Note: This method does not check for convergence of the
     iterative solution!
  */

  for (nu = 0; nu < iter_steps; nu++)
    {
      int s = GSL_ODEIV_FN_EVAL(sys, t + h, y, k);

      if (s != GSL_SUCCESS)
	{
	  return s;
	}

      for (i=0; i<dim; i++)
	{
	  y[i] = y0[i] + h * k[i];
	}
    }

  return GSL_SUCCESS;
}

static int
gear1_apply(void * vstate,
            size_t dim,
            double t,
            double h,
            double y[],
            double yerr[],
            const double dydt_in[],
            double dydt_out[],
            const gsl_odeiv_system * sys)
{
  gear1_state_t *state = (gear1_state_t *) vstate;

  size_t i;

  double *y0 = state->y0;
  double *y0_orig = state->y0_orig;
  double *y_onestep = state->y_onestep;

  /* initialization */
  DBL_MEMCPY(y0, y, dim);

  /* Save initial values for possible failures */
  DBL_MEMCPY (y0_orig, y, dim);

  /* First traverse h with one step (save to y_onestep) */
  DBL_MEMCPY (y_onestep, y, dim);

  {
    int s = gear1_step (y_onestep, state, h, t, dim, sys);

    if (s != GSL_SUCCESS)
      {
        return s;
      }
  }

 /* Then with two steps with half step length (save to y) */
  {
    int s = gear1_step (y, state, h / 2.0, t, dim, sys);

    if (s != GSL_SUCCESS)
      {
        /* Restore original y vector */
        DBL_MEMCPY (y, y0_orig, dim);
        return s;
      }
  }

  DBL_MEMCPY (y0, y, dim);

  {
    int s = gear1_step (y, state, h / 2.0, t + h / 2.0, dim, sys);

    if (s != GSL_SUCCESS)
      {
        /* Restore original y vector */
        DBL_MEMCPY (y, y0_orig, dim);
        return s;
      }
  }

  /* Cleanup update */

  if (dydt_out != NULL)
    {
      int s = GSL_ODEIV_FN_EVAL (sys, t + h, y, dydt_out);

      if (s != GSL_SUCCESS)
        {
          /* Restore original y vector */
          DBL_MEMCPY (y, y0_orig, dim);
          return s;
        }
    }

  /* Error estimation */

  for (i = 0; i < dim; i++)
    {
      yerr[i] = 4.0 * (y[i] - y_onestep[i]);
    }

  return GSL_SUCCESS;
}

static int
gear1_reset (void *vstate, size_t dim)
{
  gear1_state_t *state = (gear1_state_t *) vstate;

  DBL_ZERO_MEMSET (state->y_onestep, dim);
  DBL_ZERO_MEMSET (state->y0_orig, dim);
  DBL_ZERO_MEMSET (state->y0, dim);
  DBL_ZERO_MEMSET (state->k, dim);
  return GSL_SUCCESS;
}

static unsigned int
gear1_order (void *vstate)
{
  gear1_state_t *state = (gear1_state_t *) vstate;
  state = 0; /* prevent warnings about unused parameters */
  return 1;
}

static void
gear1_free (void *vstate)
{
  gear1_state_t *state = (gear1_state_t *) vstate;
  free (state->y_onestep);
  free (state->y0_orig);
  free (state->y0);
  free (state->k);
  free (state);
}

static const gsl_odeiv_step_type gear1_type = { "gear1",        /* name */
  0,                            /* can use dydt_in? */
  1,                            /* gives exact dydt_out? */
  &gear1_alloc,
  &gear1_apply,
  &gear1_reset,
  &gear1_order,
  &gear1_free
};

const gsl_odeiv_step_type *gsl_odeiv_step_gear1 = &gear1_type;

ext: c

max_stars_repo_name: ielomariala/Hex-Game

max_issues_repo_name: Brian-ning/HMNE

max_forks_repo_name: Brian-ning/HMNE


text:
/*
** generate permutation table using exchangeability information
**
** G.Lohmann, July 2018
*/
#include <viaio/Vlib.h>
#include <viaio/file.h>
#include <viaio/option.h>
#include <viaio/mu.h>

#include <gsl/gsl_matrix.h>
#include <gsl/gsl_vector.h>
#include <gsl/gsl_cdf.h>
#include <gsl/gsl_errno.h>
#include <gsl/gsl_math.h>
#include <gsl/gsl_rng.h>
#include <gsl/gsl_randist.h>

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

extern gsl_matrix *XRead2ndLevel(VString);
extern void XReadExchange(VString filename,int *exchange,int n);
extern gsl_matrix *PseudoInv(gsl_matrix *A,gsl_matrix *B);

void genperm(long seed,int *exchange,int signswitch,
	     int **permtable,int **signtable,int nimages,int numperm,gsl_vector *contrast)
{
  int i,j,k,m,kmax;

  /* ini permutation and sign-switching tables */
  for (i = 0; i < numperm; i++) {
    permtable[i] = (int *) VCalloc(nimages,sizeof(int));
    signtable[i] = (int *) VCalloc(nimages,sizeof(int));
    for (j=0; j<nimages; j++) permtable[i][j] = j;
    for (j=0; j<nimages; j++) signtable[i][j] = 1;
  }

  /* ini random number generator */
  gsl_rng_env_setup();
  const gsl_rng_type *T = gsl_rng_default;
  gsl_rng *rx = gsl_rng_alloc(T);
  gsl_rng_set(rx,(unsigned long int)seed);

  /* sign switching only */
  if (signswitch >= 0) {
    for (i=0; i<numperm; i++) {
      for (j=0; j<nimages; j++) {
	if (gsl_ran_bernoulli (rx,(double)0.5) == 1) signtable[i][j] = -1;
      }
    }
    return;  /* stop here, no permutations */
  }

  /* get max number of exchange groups */
  kmax=0;
  for (j=0; j<nimages; j++) {
    if (exchange[j] > kmax) kmax = exchange[j];
  }

  /* permutations of rows of the design matrix */
  int *base = (int *) VCalloc(nimages,sizeof(int));
  for (i = 0; i < numperm; i++) {

    /* for each group */
    for (k=0; k<=kmax; k++) {
      m = 0;
      for (j=0; j<nimages; j++) {
	if (exchange[j] == k) {
	  base[m] = j;
	  m++;
	}
      }
      if (m < 1) continue;

      /* random shuffle within group */
      gsl_ran_shuffle (rx, base, m, sizeof (int));

      m = 0;
      for (j=0; j<nimages; j++) {
	if (exchange[j] == k) {
	  permtable[i][j] = base[m];
	  m++;
	}
      }
    }
  }
  VFree(base);
}

/* for zero variance replace permutations with sign switching (only for non-nuisance columns) */
int SignSwitch (gsl_matrix *X,gsl_vector *contrast,int *permflag)
{
  int i,j;
  double s1=0,s2=0,nx=0,u=0,mean=0,var=0;

  int signswitch = -1;
  for (j=0; j<X->size2; j++) {
    if (permflag[j] == 0) continue;
    if (fabs(contrast->data[j]) < TINY) continue;
    s1 = s2 = nx = 0;
    for (i=0; i<X->size1; i++) {
      u = gsl_matrix_get(X,i,j);
      s1 += u;
      s2 += u*u;
      nx++;
    }
    mean = s1/nx;
    var = (s2 - nx * mean * mean) / (nx - 1.0);
    if (var < TINY && signswitch >= 0)
      VError(" Implausible design, two columns have zero variance: %d %d",signswitch+1,j+1);
    if (var < TINY) signswitch = j;
  }
  return signswitch;
}

max_stars_repo_name: zrajna/lipsia

max_issues_repo_name: zrajna/lipsia

max_forks_repo_name: zrajna/lipsia


text:
#ifndef _COSMOLOGY_H_
#define _COSMOLOGY_H_

#include <gsl/gsl_integration.h>

class cosmology{
    double Om_M, Om_L, Om_b, Om_c, tau, T_CMB, h;

    double E(double z);

    static double E_inv(double z, void *params);

    static double rd_int(double z, void *params);

    static double rz(double z, void *params);

    double Theta();

    double k_eq();

    public:
        cosmology(double H_0 = 70.0, double OmegaM = 0.3, double OmegaL = 0.7, double Omegab = 0.04,
                  double Omegac = 0.26, double Tau = 0.066, double TCMB = 2.718);

        double Omega_M();

        double Omega_L();

        double Omega_bh2();

        double Omega_ch2();

        double h_param();

        double H0();

        double H(double z);

        double D_A(double z);

        double D_V(double z);

        double z_eq();

        double z_d();

        double R(double z);

        double r_d();

        double comoving_distance(double z, gsl_integration_workspace *w);

};

#endif


text:
/* sys/infnan.c
 *
 * Copyright (C) 2001, 2004, 2007, 2010 Brian Gough
 *
 * 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 <math.h>

#if HAVE_IEEEFP_H
#include <ieeefp.h>
#endif

#include <gsl/gsl_sys.h>

double gsl_nan (void)
{
  return gsl_fdiv (0.0, 0.0);
}

double gsl_posinf (void)
{
  return gsl_fdiv (+1.0, 0.0);
}

double gsl_neginf (void)
{
  return gsl_fdiv (-1.0, 0.0);
}

int gsl_isnan (const double x);
int gsl_isinf (const double x);
int gsl_finite (const double x);

#if defined(_MSC_VER) /* Microsoft Visual C++ */
#include <float.h>
int
gsl_isnan (const double x)
{
  return _isnan(x);
}

int
gsl_isinf (const double x)
{
  int fpc = _fpclass(x);

  if (fpc == _FPCLASS_PINF)
    return +1;
  else if (fpc == _FPCLASS_NINF)
    return -1;
  else
    return 0;
}

int
gsl_finite (const double x)
{
  return _finite(x);
}
#else

# if HAVE_DECL_ISFINITE
int
gsl_finite (const double x)
{
  return isfinite(x);
}
# elif HAVE_DECL_FINITE
int
gsl_finite (const double x)
{
  return finite(x);
}
# elif HAVE_IEEE_COMPARISONS
int
gsl_finite (const double x)
{
  const double y = x - x;
  int status = (y == y);
  return status;
}
# else
# error "cannot define gsl_finite without HAVE_DECL_FINITE or HAVE_IEEE_COMPARISONS"
# endif

# if HAVE_DECL_ISNAN
int
gsl_isnan (const double x)
{
  return isnan(x);
}
#elif HAVE_IEEE_COMPARISONS
int
gsl_isnan (const double x)
{
  int status = (x != x);
  return status;
}
# else
# error "cannot define gsl_isnan without HAVE_DECL_ISNAN or HAVE_IEEE_COMPARISONS"
# endif

# if HAVE_DECL_ISINF
int
gsl_isinf (const double x)
{
  /* isinf(3): In glibc 2.01 and earlier, isinf() returns a
     non-zero value (actually: 1) if x is an infinity (positive or
     negative).  (This is all that C99 requires.) */

  if (isinf(x))
    {
      return (x > 0) ? 1 : -1;
    }
  else
    {
      return 0;
    }
}
# else

int
gsl_isinf (const double x)
{
  if (! gsl_finite(x) && ! gsl_isnan(x))
    {
      return (x > 0 ? +1 : -1);
    }
  else
    {
      return 0;
    }
}

# endif
#endif

max_stars_repo_licenses: Net-SNMP Xnet

max_issues_repo_licenses: Net-SNMP Xnet

max_forks_repo_licenses: Net-SNMP Xnet


text:
/**
 *
 * @file zgesdd.c
 *
 *  PLASMA computational routines
 *  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
 * @date 2010-11-15
 * @precisions normal z -> s d c
 *
 **/
#include <lapacke.h>
#include "common.h"

#define COMPLEX
/***************************************************************************//**
 *
 * @ingroup PLASMA_Complex64_t
 *
 *  PLASMA_zgesdd - computes the singular value decomposition (SVD) of a complex
 *  M-by-N matrix A, optionally computing the left and/or right singular
 *  vectors, by using divide-and-conquer method. The SVD is written
 *
 *       A = U * SIGMA * conjugate-transpose(V)
 *
 *  where SIGMA is an M-by-N matrix which is zero except for its
 *  min(m,n) diagonal elements, U is an M-by-M unitary matrix, and
 *  V is an N-by-N unitary matrix.  The diagonal elements of SIGMA
 *  are the singular values of A; they are real and non-negative, and
 *  are returned in descending order.  The first min(m,n) columns of
 *  U and V are the left and right singular vectors of A.
 *
 *  Note that the routine returns VT = V**H, not V.
 *
 *  The divide and conquer algorithm makes very mild assumptions about
 *  floating point arithmetic. It will work on machines with a guard
 *  digit in add/subtract, or on those binary machines without guard
 *  digits which subtract like the Cray X-MP, Cray Y-MP, Cray C-90, or
 *  Cray-2. It could conceivably fail on hexadecimal or decimal machines
 *  without guard digits, but we know of none.
 *
 *  NOTE that the input parameter of this function differ from the
 *  LAPACK ZGESDD. We have jobu and jobvt instead of only jobz first
 *  to be consistent with ZGESVD and second to allow a computation
 *  of only U or VT.
 *******************************************************************************
 *
 * @param[in] jobu
 *          Specifies options for computing all or part of the matrix U.
 *          Intended usage:
 *          = PlasmaVec   = 'A'(lapack):  all M columns of U are returned
 *                        in array U;
 *          = PlasmaNoVec = 'N':  no columns of U (no left singular vectors)
 *                        are computed.
 *          = PlasmaSVec  = 'S': the first min(m,n) columns of U (the left
 *                        singular vectors) are returned in the array U;
 *                        NOT SUPPORTTED YET
 *          = PlasmaOVec  = 'O': the first min(m,n) columns of U (the left
 *                        singular vectors) are overwritten on the array A;
 *                        NOT SUPPORTTED YET
 *
 * @param[in] jobvt
 *          Specifies options for computing all or part of the matrix V**H.
 *          Intended usage:
 *          = PlasmaVec   = 'A'(lapack): all N rows of V**H are returned
 *                        in the array VT;
 *          = PlasmaNoVec = 'N': no rows of V**H (no right singular vectors)
 *                        are computed.
 *          = PlasmaSVec  = 'S': the first min(m,n) rows of V**H (the right
 *                        singular vectors) are returned in the array VT;
 *                        NOT SUPPORTTED YET
 *          = PlasmaOVec  = 'O': the first min(m,n) rows of V**H (the right
 *                        singular vectors) are overwritten on the array A;
 *                        NOT SUPPORTTED YET
 *
 *          Note: jobu and jobvt cannot both be PlasmaOVec.
 *
 * @param[in] M
 *          The number of rows of the matrix A. M >= 0.
 *
 * @param[in] N
 *          The number of columns of the matrix A. N >= 0.
 *
 * @param[in,out] A
 *          On entry, the M-by-N matrix A.
 *          On exit,
 *          if JOBU = 'O',  A is overwritten with the first min(m,n)
 *                          columns of U (the left singular vectors,
 *                          stored columnwise);
 *          if JOBVT = 'O', A is overwritten with the first min(m,n)
 *                          rows of V**H (the right singular vectors,
 *                          stored rowwise);
 *          if JOBU .ne. 'O' and JOBVT .ne. 'O', the contents of A
 *                          are destroyed.
 *
 * @param[in] LDA
 *          The leading dimension of the array A. LDA >= max(1,M).
 *
 * @param[out] S
 *          The double precision singular values of A, sorted so that S(i) >= S(i+1).
 *
 * @param[in, out] descT
 *          On entry, descriptor as return by PLASMA_Alloc_Workspace_zgesdd
 *          On exit, contains auxiliary factorization data.
 *
 * @param[out] U
 *          (LDU,M) if JOBU = 'A' or (LDU,min(M,N)) if JOBU = 'S'.
 *          If JOBU = 'A', U contains the M-by-M unitary matrix U;
 *          if JOBU = 'S', U contains the first min(m,n) columns of U
 *          (the left singular vectors, stored columnwise);
 *          if JOBU = 'N' or 'O', U is not referenced.
 *
 * @param[in] LDU
 *          The leading dimension of the array U.  LDU >= 1; if
 *          JOBU = 'S' or 'A', LDU >= M.
 *
 * @param[out] VT
 *         If JOBVT = 'A', VT contains the N-by-N unitary matrix
 *         V**H;
 *         if JOBVT = 'S', VT contains the first min(m,n) rows of
 *         V**H (the right singular vectors, stored rowwise);
 *         if JOBVT = 'N' or 'O', VT is not referenced.
 *
 * @param[in] LDVT
 *         The leading dimension of the array VT.  LDVT >= 1; if
 *         JOBVT = 'A', LDVT >= N; if JOBVT = 'S', LDVT >= min(M,N).
 *
 *******************************************************************************
 *
 * @return
 *          \retval PLASMA_SUCCESS successful exit
 *          \retval <0 if -i, the i-th argument had an illegal value
 *
 *******************************************************************************
 *
 * @sa PLASMA_zgesdd_Tile
 * @sa PLASMA_zgesdd_Tile_Async
 * @sa PLASMA_cgesdd
 * @sa PLASMA_dgesdd
 * @sa PLASMA_sgesdd
 *
 ******************************************************************************/
int PLASMA_zgesdd(PLASMA_enum jobu, PLASMA_enum jobvt, int M, int N,
                  PLASMA_Complex64_t *A, int LDA,
                  double *S,
                  PLASMA_desc *descT,
                  PLASMA_Complex64_t *U, int LDU,
                  PLASMA_Complex64_t *VT, int LDVT)
{
    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_fatal_error("PLASMA_zgesdd", "PLASMA not initialized");
        return PLASMA_ERR_NOT_INITIALIZED;
    }

    /* Check input arguments */
    if (jobu != PlasmaNoVec  && jobu !=PlasmaVec) {
        plasma_error("PLASMA_zgesdd", "illegal value of jobu");
        return -1;
    }
    if (jobvt != PlasmaNoVec && jobvt != PlasmaVec) {
        plasma_error("PLASMA_zgesdd", "illegal value of jobvt");
        return -2;
    }
    if (jobvt != jobu) {
        plasma_error("PLASMA_zgesdd", "in this version: jobu should be equal jobvt");
        return -22;
    }
    if (M < 0) {
        plasma_error("PLASMA_zgesdd", "illegal value of M");
        return -3;
    }
    if (N < 0) {
        plasma_error("PLASMA_zgesdd", "illegal value of N");
        return -4;
    }
    if (LDA < max(1, M)) {
        plasma_error("PLASMA_zgesdd", "illegal value of LDA");
        return -6;
    }
    if (LDU < 1) {
        plasma_error("PLASMA_zgesdd", "illegal value of LDU");
        return -9;
    }
    if (LDVT < 1) {
        plasma_error("PLASMA_zgesdd", "illegal value of LDVT");
        return -11;
    }
    /* Quick return */
    if (min(M, N) == 0) {
        return PLASMA_SUCCESS;
    }

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

    /* Set MT, NT */
    NB = PLASMA_NB;

    plasma_sequence_create(plasma, &sequence);

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

    /* Call the tile interface */
    PLASMA_zgesdd_Tile_Async(jobu, jobvt, &descA, S, descT, U, LDU, VT, LDVT, 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_zgesdd_Tile - computes the singular value decomposition (SVD) of a complex
 *  M-by-N matrix A, optionally computing the left and/or right singular
 *  vectors.
 *  Tile equivalent of PLASMA_zgesdd().
 *  Operates on matrices stored by tiles.
 *  All matrices are passed through descriptors.
 *  All dimensions are taken from the descriptors.
 *
 *******************************************************************************
 *
 * @param[in] jobu
 *          Specifies options for computing all or part of the matrix U.
 *          Intended usage:
 *          = PlasmaVec   = 'A'(lapack):  all M columns of U are returned
 *                        in array U;
 *          = PlasmaNoVec = 'N':  no columns of U (no left singular vectors)
 *                        are computed.
 *          = PlasmaSVec  = 'S': the first min(m,n) columns of U (the left
 *                        singular vectors) are returned in the array U;
 *                        NOT SUPPORTTED YET
 *          = PlasmaOVec  = 'O': the first min(m,n) columns of U (the left
 *                        singular vectors) are overwritten on the array A;
 *                        NOT SUPPORTTED YET
 *
 * @param[in] jobvt
 *          Specifies options for computing all or part of the matrix V**H.
 *          Intended usage:
 *          = PlasmaVec   = 'A'(lapack): all N rows of V**H are returned
 *                        in the array VT;
 *          = PlasmaNoVec = 'N': no rows of V**H (no right singular vectors)
 *                        are computed.
 *          = PlasmaSVec  = 'S': the first min(m,n) rows of V**H (the right
 *                        singular vectors) are returned in the array VT;
 *                        NOT SUPPORTTED YET
 *          = PlasmaOVec  = 'O': the first min(m,n) rows of V**H (the right
 *                        singular vectors) are overwritten on the array A;
 *                        NOT SUPPORTTED YET
 *
 *          Note: jobu and jobvt cannot both be PlasmaOVec.
 *
 * @param[in,out] A
 *          On entry, the M-by-N matrix A.
 *          On exit,
 *          if JOBU = 'O',  A is overwritten with the first min(m,n)
 *                          columns of U (the left singular vectors,
 *                          stored columnwise);
 *          if JOBVT = 'O', A is overwritten with the first min(m,n)
 *                          rows of V**H (the right singular vectors,
 *                          stored rowwise);
 *          if JOBU .ne. 'O' and JOBVT .ne. 'O', the contents of A
 *                          are destroyed.
 *
 * @param[out] S
 *          The singular values of A, sorted so that S(i) >= S(i+1).
 *
 * @param[in, out] T
 *          On entry, descriptor as return by PLASMA_Alloc_Workspace_zgesdd
 *          On exit, contains auxiliary factorization data.
 *
 * @param[out] U
 *          (LDU,M) if JOBU = 'A' or (LDU,min(M,N)) if JOBU = 'S'.
 *          If JOBU = 'A', U contains the M-by-M unitary matrix U;
 *          if JOBU = 'S', U contains the first min(m,n) columns of U
 *          (the left singular vectors, stored columnwise);
 *          if JOBU = 'N' or 'O', U is not referenced.
 *
 * @param[in] LDU
 *          The leading dimension of the array U.  LDU >= 1; if
 *          JOBU = 'S' or 'A', LDU >= M.
 *
 * @param[out] VT
 *         If JOBVT = 'A', VT contains the N-by-N unitary matrix
 *         V**H;
 *         if JOBVT = 'S', VT contains the first min(m,n) rows of
 *         V**H (the right singular vectors, stored rowwise);
 *         if JOBVT = 'N' or 'O', VT is not referenced.
 *
 * @param[in] LDVT
 *         The leading dimension of the array VT.  LDVT >= 1; if
 *         JOBVT = 'A', LDVT >= N; if JOBVT = 'S', LDVT >= min(M,N).
 *
 *******************************************************************************
 *
 * @return
 *          \return PLASMA_SUCCESS successful exit
 *
 *******************************************************************************
 *
 * @sa PLASMA_zgesdd
 * @sa PLASMA_zgesdd_Tile_Async
 * @sa PLASMA_cgesdd_Tile
 * @sa PLASMA_dgesdd_Tile
 * @sa PLASMA_sgesdd_Tile
 *
 ******************************************************************************/
int PLASMA_zgesdd_Tile(PLASMA_enum jobu, PLASMA_enum jobvt,
                       PLASMA_desc *A,
                       double *S,
                       PLASMA_desc *T,
                       PLASMA_Complex64_t *U, int LDU,
                       PLASMA_Complex64_t *VT, int LDVT)
{
    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_zgesdd_Tile", "PLASMA not initialized");
        return PLASMA_ERR_NOT_INITIALIZED;
    }
    plasma_sequence_create(plasma, &sequence);
    PLASMA_zgesdd_Tile_Async(jobu, jobvt, A, S, T, U, LDU, VT, LDVT, sequence, &request);
    plasma_dynamic_sync();
    status = sequence->status;
    plasma_sequence_destroy(plasma, sequence);
    return status;
}

/***************************************************************************//**
 *
 * @ingroup PLASMA_Complex64_t_Tile_Async
 *
 *  PLASMA_zgesdd_Tile_Async - computes the singular value decomposition (SVD) of a complex
 *  M-by-N matrix A, optionally computing the left and/or right singular
 *  vectors.
 *  Non-blocking equivalent of PLASMA_zgesdd_Tile().
 *  May return before the computation is finished.
 *  Allows for pipelining of operations at runtime.
 *
 *******************************************************************************
 *
 * @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_zgesdd
 * @sa PLASMA_zgesdd_Tile
 * @sa PLASMA_cgesdd_Tile_Async
 * @sa PLASMA_dgesdd_Tile_Async
 * @sa PLASMA_sgesdd_Tile_Async
 *
 ******************************************************************************/
int PLASMA_zgesdd_Tile_Async(PLASMA_enum jobu, PLASMA_enum jobvt,
                             PLASMA_desc *A,
                             double *S,
                             PLASMA_desc *T,
                             PLASMA_Complex64_t *U, int LDU,
                             PLASMA_Complex64_t *VT, int LDVT,
                             PLASMA_sequence *sequence, PLASMA_request *request)
{
    PLASMA_desc descA;
    PLASMA_desc descT;
    PLASMA_desc descU, descVT;
    PLASMA_Complex64_t *AB;
    double *E;
    int M;
    int N;
    int MINMN;
    int NB;
    int LDAB;
    int i;
    int status;

    plasma_context_t *plasma;
    plasma = plasma_context_self();

    if (plasma == NULL) {
        plasma_fatal_error("PLASMA_zgesdd_Tile_Async", "PLASMA not initialized");
        return PLASMA_ERR_NOT_INITIALIZED;
    }
    if (sequence == NULL) {
        plasma_fatal_error("PLASMA_zgesdd_Tile_Async", "NULL sequence");
        return PLASMA_ERR_UNALLOCATED;
    }
    if (request == NULL) {
        plasma_fatal_error("PLASMA_zgesdd_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_zgesdd_Tile_Async", "invalid first descriptor");
        return plasma_request_fail(sequence, request, PLASMA_ERR_ILLEGAL_VALUE);
    } else {
        descA = *A;
    }
    if (plasma_desc_check(T) != PLASMA_SUCCESS) {
        plasma_error("PLASMA_zgesdd_Tile_Async", "invalid fourth descriptor");
        return plasma_request_fail(sequence, request, PLASMA_ERR_ILLEGAL_VALUE);
    } else {
        descT = *T;
    }
    /* Check input arguments */
    if (jobu != PlasmaNoVec  && jobu != PlasmaVec) {
        plasma_error("PLASMA_zgesdd_Tile_Async", "illegal value of jobu");
        return PLASMA_ERR_NOT_SUPPORTED;
    }
    if (jobvt != PlasmaNoVec && jobvt != PlasmaVec) {
        plasma_error("PLASMA_zgesdd_Tile_Async", "illegal value of jobvt");
        return PLASMA_ERR_NOT_SUPPORTED;
    }
    if (jobvt != jobu) {
        plasma_error("PLASMA_zgesdd_Tile_Async", "in this version: jobu should be equal jobvt");
        return PLASMA_ERR_NOT_SUPPORTED;
    }
    if (descA.nb != descA.mb) {
        plasma_error("PLASMA_zgesdd_Tile_Async", "only square tiles supported");
        return plasma_request_fail(sequence, request, PLASMA_ERR_ILLEGAL_VALUE);
    }

    #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();
    #endif
    PLASMA_enum uplo = descA.m >= descA.n ? PlasmaUpper : PlasmaLower;
    M     = descA.m;
    N     = descA.n;
    MINMN = min(M,N);
    NB    = min(descA.mb,MINMN);
    LDAB  = 3*NB+1;
    /*=======================================
     *  case M<NB or N<NB call lapack
     *=======================================*/
    if( ( M<= NB) || ( N <=NB ) ){
        /* convert the tile descA to lapack A */
        if ( PLASMA_TRANSLATION == PLASMA_OUTOFPLACE ) {
            AB = (PLASMA_Complex64_t *) plasma_shared_alloc(plasma, M*N, PlasmaComplexDouble);
            if (AB == NULL) {
                plasma_error("PLASMA_zgesvd_Tile_Async", "plasma_shared_alloc(AB-0-) failed");
                plasma_shared_free(plasma, AB);
                return PLASMA_ERR_OUT_OF_RESOURCES;
            }
            plasma_zooptile2lap( descA, AB, NB, NB, M, N,  sequence, request);
        } else {
            AB = descA.mat;
            plasma_ziptile2lap( descA, AB, NB, NB, M, N,  sequence, request);
        }
        plasma_dynamic_sync();
        /*=======================================
         *  calling LAPACK ZGESVD
         *=======================================*/
        #if defined(ENABLE_TIMER)
        plasma_dynamic_sync();
        timelpk   = PLASMA_Wtime();
        #endif
        plasma_setlapack_multithreads(plasma->world_size);
        /* call SVD solver using lapack routine  */
        PLASMA_enum jobz  = jobu  == PlasmaVec ? PlasmaAllVec : PlasmaNoVec;
        status = LAPACKE_zgesdd(LAPACK_COL_MAJOR, lapack_const(jobz),
                                M, N, AB, M, S, U, LDU, VT, LDVT);
        if(status != 0){
            plasma_error("PLASMA_zgesvd","ZGESVD");
        }
        sequence->status = status;
        plasma_setlapack_sequential(plasma);
        #if defined(ENABLE_TIMER)
        timelpk = PLASMA_Wtime()-timelpk;
        printf("  Finish Eigensolver-lpkonly timing= %lf  threads %d\n" ,timelpk, plasma->world_size);
        #endif
        /*=======================================
         *  END of calling SVD solver
         *=======================================*/
        /* convert the lapack to the tile descA */
        if ( PLASMA_TRANSLATION == PLASMA_OUTOFPLACE ) {
            //plasma_zooplap2tile_noalloc( descA, AB, NB, NB,  M, N, 0, 0, M, N, sequence, request);
            plasma_parallel_call_5( plasma_pzlapack_to_tile,
            PLASMA_Complex64_t*, AB,
            int,                 M,
            PLASMA_desc,         descA,
            PLASMA_sequence*,    sequence,
            PLASMA_request*,     request);
            plasma_dynamic_sync();
            free(AB);
        } else {
            plasma_ziplap2tile( descA, AB, NB, NB,  M, N, 0, 0, M, N,
                                sequence, request);
        }
        plasma_dynamic_sync();
        return PLASMA_SUCCESS;
    }
    /*=======================================
     *  END OF case M<NB or N<NB
     *=======================================*/
    /*
     * Allocate workspace for band storage of the band matrix A
     * AB looks like:
     *       __________________________________
     * NB   |               zero               |
     *       ----------------------------------
     * NB+1 |               band A             |
     *       ----------------------------------
     * NB   |_______________zero_______________|
     *
     * */
    AB = (PLASMA_Complex64_t *)plasma_shared_alloc(plasma, LDAB*MINMN, PlasmaComplexDouble);
    memset( AB, 0, LDAB * MINMN * sizeof(PLASMA_Complex64_t) );
    if (AB == NULL) {
        plasma_error("PLASMA_zgesdd_Tile_Async", "plasma_shared_alloc(AB) failed");
        plasma_shared_free(plasma, AB);
        return PLASMA_ERR_OUT_OF_RESOURCES;
    }
    E = (double *)plasma_shared_alloc(plasma, MINMN, PlasmaRealDouble);
    if (E == NULL) {
        plasma_error("PLASMA_zgesdd_Tile_Async", "plasma_shared_alloc(E) failed");
        plasma_shared_free(plasma, E);
        return PLASMA_ERR_OUT_OF_RESOURCES;
    }
    /*=======================================
     *  calling Reduction from DENSE to BAND
     *  then convert matrix to band form
     *=======================================*/
    //plasma_dynamic_sync();
    #if defined(ENABLE_TIMER)
    timeB   = PLASMA_Wtime();
    #endif
    /*
     * Reduction to BAND bidiagonal form
     * May be further optimized using the algo described in Trefethen
     */
    /* if (plasma->householder == PLASMA_FLAT_HOUSEHOLDER) { */
        plasma_dynamic_call_4(plasma_pzgebrd_ge2gb,
            PLASMA_desc, descA,
            PLASMA_desc, descT,
            PLASMA_sequence*, sequence,
            PLASMA_request*, request);
    /* } */
    /* else { */
    /*     plasma_dynamic_call_4(plasma_pzgebrd_ge2gb_rh, */
    /*         PLASMA_desc, descA, */
    /*         PLASMA_desc, descT, */
    /*         PLASMA_sequence*, sequence, */
    /*         PLASMA_request*, request); */
    /* } */
    //plasma_dynamic_sync();
    plasma_dynamic_call_6( plasma_pzgbcpy_t2bl,
        PLASMA_enum, descA.m >= descA.n ? PlasmaUpper : PlasmaLower,
        PLASMA_desc, descA,
        PLASMA_Complex64_t*, &(AB[NB]),
        int, LDAB,
        PLASMA_sequence*, sequence,
        PLASMA_request*, request);
    plasma_dynamic_sync();
    status = sequence->status;
    if (status != PLASMA_SUCCESS) {
        plasma_error("PLASMA_zgesdd","pzgebrd_ge2gb+pzcopy");
        return status;
    }
    #if defined(ENABLE_TIMER)
    timeB   = PLASMA_Wtime()-timeB;
    printf("\n  Finish Band red    timing= %lf \n",timeB);
    #endif
    /*=======================================
     *  END of calling Reduction to BAND
     *=======================================*/
    /*=======================================
     *  calling Reduction from BAND to bidiag
     *=======================================*/
    PLASMA_Complex64_t *VQ2   = NULL;
    PLASMA_Complex64_t *VP2   = NULL;
    PLASMA_Complex64_t *TAUQ2 = NULL;
    PLASMA_Complex64_t *TAUP2 = NULL;
    PLASMA_Complex64_t *TQ2   = NULL;
    PLASMA_Complex64_t *TP2   = NULL;
    int Vblksiz, blkcnt, LDT, LDV;
    int WANTZ   = 0;
    /* compute fraction of singular vectors */
    /* for future use when a portion of the eigenvectors are requested */
    int ME      = M;
    int NE      = N;
    int MINMNE  = MINMN;
    /*
    int fraction = 100;
    if((fraction<100)&&(fraction>0)) {
        MINMNE  = plasma_ceildiv(fraction*MINMN,100);
        ME      = MINMNE;
        NE      = MINMNE;
    }
    */
    if( jobu == PlasmaNoVec )
        WANTZ=0;
    else
        WANTZ=2;
    /* Vblksiz correspond to the blocking used when applying V2 to the matrix U
     * it is similar to IB in LAPACK ZUNMQR.
     * blkcnt is the number of diamond 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 diamond in a tile storage meaning
     * that V2_1 is stored then V2_2,..., V2_blkcnt.
     * blkcnt is the number of diamond.
     * */
    Vblksiz = min(NB,48);
    LDT     = Vblksiz;
    /* data for U */
    if( jobu == PlasmaVec ) {
        findVTsiz(MINMN, NB, Vblksiz, &blkcnt, &LDV);
        TAUQ2   = (PLASMA_Complex64_t *) plasma_shared_alloc(plasma,     blkcnt*Vblksiz, PlasmaComplexDouble);
        VQ2     = (PLASMA_Complex64_t *) plasma_shared_alloc(plasma, LDV*blkcnt*Vblksiz, PlasmaComplexDouble);
        TQ2     = (PLASMA_Complex64_t *) plasma_shared_alloc(plasma, LDT*blkcnt*Vblksiz, PlasmaComplexDouble);
        if ( (TAUQ2 == NULL) || (VQ2 == NULL) || (TQ2 == NULL) ) {
            plasma_error("PLASMA_zgesdd", "plasma_shared_alloc() failed");
            plasma_shared_free(plasma, TAUQ2);
            plasma_shared_free(plasma, VQ2);
            plasma_shared_free(plasma, TQ2);
            return PLASMA_ERR_OUT_OF_RESOURCES;
        }
        memset(TAUQ2, 0,     blkcnt*Vblksiz*sizeof(PLASMA_Complex64_t));
        memset(VQ2,   0, LDV*blkcnt*Vblksiz*sizeof(PLASMA_Complex64_t));
        memset(TQ2,   0, LDT*blkcnt*Vblksiz*sizeof(PLASMA_Complex64_t));
    }
    else {
        TAUQ2   = (PLASMA_Complex64_t *) plasma_shared_alloc(plasma, 2*MINMN, PlasmaComplexDouble);
        VQ2     = (PLASMA_Complex64_t *) plasma_shared_alloc(plasma, 2*MINMN, PlasmaComplexDouble);
        if ( (TAUQ2 == NULL) || (VQ2 == NULL) ) {
            plasma_error("PLASMA_zgesdd", "plasma_shared_alloc() failed");
            plasma_shared_free(plasma, TAUQ2);
            plasma_shared_free(plasma, VQ2);
            return PLASMA_ERR_OUT_OF_RESOURCES;
        }
        memset(TAUQ2, 0, 2*MINMN*sizeof(PLASMA_Complex64_t));
        memset(VQ2,   0, 2*MINMN*sizeof(PLASMA_Complex64_t));
    }
    /* data for VT */
    if( jobvt == PlasmaVec ) {
        findVTsiz(MINMN, NB, Vblksiz, &blkcnt, &LDV);
        TAUP2   = (PLASMA_Complex64_t *) plasma_shared_alloc(plasma,     blkcnt*Vblksiz, PlasmaComplexDouble);
        VP2     = (PLASMA_Complex64_t *) plasma_shared_alloc(plasma, LDV*blkcnt*Vblksiz, PlasmaComplexDouble);
        TP2     = (PLASMA_Complex64_t *) plasma_shared_alloc(plasma, LDT*blkcnt*Vblksiz, PlasmaComplexDouble);
        if ( (TAUP2 == NULL) || (VP2 == NULL) || (TP2 == NULL) ) {
            plasma_error("PLASMA_zgesdd", "plasma_shared_alloc() failed");
            plasma_shared_free(plasma, TAUP2);
            plasma_shared_free(plasma, VP2);
            plasma_shared_free(plasma, TP2);
            return PLASMA_ERR_OUT_OF_RESOURCES;
        }
        memset(TAUP2, 0,     blkcnt*Vblksiz*sizeof(PLASMA_Complex64_t));
        memset(VP2,   0, LDV*blkcnt*Vblksiz*sizeof(PLASMA_Complex64_t));
        memset(TP2,   0, LDT*blkcnt*Vblksiz*sizeof(PLASMA_Complex64_t));
    }
    else {
        TAUP2   = (PLASMA_Complex64_t *) plasma_shared_alloc(plasma, 2*MINMN, PlasmaComplexDouble);
        VP2     = (PLASMA_Complex64_t *) plasma_shared_alloc(plasma, 2*MINMN, PlasmaComplexDouble);
        if ( (TAUP2 == NULL) || (VP2 == NULL) ) {
            plasma_error("PLASMA_zgesdd", "plasma_shared_alloc() failed");
            plasma_shared_free(plasma, TAUQ2);
            plasma_shared_free(plasma, VQ2);
            return PLASMA_ERR_OUT_OF_RESOURCES;
        }
        memset(TAUP2, 0, 2*MINMN*sizeof(PLASMA_Complex64_t));
        memset(VP2,   0, 2*MINMN*sizeof(PLASMA_Complex64_t));
    }
    /*=======================================
     *  calling bulge chasing
     *=======================================*/
    #if defined(ENABLE_TIMER)
    timeblg  = PLASMA_Wtime();
    #endif
    plasma_parallel_call_16(plasma_pzgebrd_gb2bd_v1,
        PLASMA_enum,         uplo,
        int,                 MINMN,
        int,                 NB,
        int,                 Vblksiz,
        PLASMA_Complex64_t*, AB,
        int,                 LDAB,
        PLASMA_Complex64_t*, VQ2,
        PLASMA_Complex64_t*, TAUQ2,
        PLASMA_Complex64_t*, VP2,
        PLASMA_Complex64_t*, TAUP2,
        double*,             S,
        double*,             E,
        int,                 WANTZ,
        int,                 WANTZ,
        PLASMA_sequence*,    sequence,
        PLASMA_request*,     request);
    /* WARNING: If plasma_pzgebrd_gb2bd is implemented through a dynamic call, don't
     * forget to synchronize */
    plasma_dynamic_sync();
    #if defined(ENABLE_TIMER)
    timeblg   = PLASMA_Wtime()-timeblg;
    printf("  Finish Bulge       timing= %lf \n" ,timeblg);
    #endif
    /*=======================================
     *  END of calling bulge chasing
     *=======================================*/
    /*=======================================
     *  calling SVD solver
     *=======================================*/
    plasma_setlapack_multithreads(plasma->world_size);
    /* call SVD solver using lapack routine for our resulting bidiag [S E] */
    /* call eigensolver using lapack routine for our resulting tridiag [D E] */
    if((jobu == PlasmaNoVec)&&(jobvt == PlasmaNoVec)){
        #if defined(ENABLE_TIMER)
        timelpk   = PLASMA_Wtime();
        #endif
        status = LAPACKE_dbdsdc(LAPACK_COL_MAJOR, lapack_const(uplo), lapack_const(PlasmaNoVec),
                                MINMN, S, E, NULL, LDU, NULL, LDVT, NULL, NULL);
        #if defined(ENABLE_TIMER)
        timelpk = PLASMA_Wtime()-timelpk;
        #endif
    } else {
        // set U and VT to zero
        memset(VT,   0, N*LDVT*sizeof(PLASMA_Complex64_t));
        memset(U,   0, M*LDU*sizeof(PLASMA_Complex64_t));
        /* Initialize U(N+1:M,N+1:M) and VT(M+1:N,M+1:N) to Identity */

        for(i=N; i<M; i++){
            U[i+i*LDU] = 1.0;
        }
         for(i=M; i<N; i++){
            VT[i+i*LDVT] = 1.0;
        }
#if defined COMPLEX
        int j;
        double *RU    = NULL;
        double *RVT   = NULL;
        RU   = (double *) plasma_shared_alloc(plasma,     MINMN*MINMN, PlasmaRealDouble);
        RVT  = (double *) plasma_shared_alloc(plasma,     MINMN*MINMN, PlasmaRealDouble);
        #if defined(ENABLE_TIMER)
        timelpk   = PLASMA_Wtime();
        #endif
        status = LAPACKE_dbdsdc(LAPACK_COL_MAJOR, lapack_const(uplo), lapack_const(PlasmaIvec),
                                MINMN, S, E, RU, MINMN, RVT, MINMN, NULL, NULL);
        #if defined(ENABLE_TIMER)
        timelpk = PLASMA_Wtime()-timelpk;
        #endif
        // copy the real U and VT matrix to the complex
        //LAPACKE_zlacp2(LAPACK_COL_MAJOR, lapack_const(uplo), MINMN, MINMN, RU, MINMN, U, LDU );
        //LAPACKE_zlacp2(LAPACK_COL_MAJOR, lapack_const(uplo), MINMN, MINMN, RVT, MINMN, VT, LDVT );
        for(j=0; j<MINMN; j++){
            for(i=0; i<MINMN; i++){
                U[i+j*LDU] = (PLASMA_Complex64_t) RU[i+j*MINMN];
            }
        }
        for(j=0; j<MINMN; j++){
            for(i=0; i<MINMN; i++){
                VT[i+j*LDVT] = (PLASMA_Complex64_t) RVT[i+j*MINMN];
            }
        }
        plasma_shared_free(plasma, RU);
        plasma_shared_free(plasma, RVT);
#else
        #if defined(ENABLE_TIMER)
        timelpk   = PLASMA_Wtime();
        #endif
        status = LAPACKE_dbdsdc(LAPACK_COL_MAJOR, lapack_const(uplo), lapack_const(PlasmaIvec),
                                MINMN, S, E, U, LDU, VT, LDVT, NULL, NULL);
        #if defined(ENABLE_TIMER)
        timelpk = PLASMA_Wtime()-timelpk;
        #endif
#endif
    }
    if(status != 0){
        plasma_error("PLASMA_zbdedc","ZBDEDC");
        /*return status;*/
    }
    sequence->status = status;
    plasma_setlapack_sequential(plasma);
    #if defined(ENABLE_TIMER)
    int fraction = (MINMNE*100)/N;
    printf("  Finish Eigensolver timing= %lf  WANTZ %d threads %d fraction %d  MINMNE %d \n" ,timelpk, WANTZ,  plasma->world_size, fraction, MINMNE);
    #endif
    /*=======================================
     *  END of calling SVD solver
     *=======================================*/
    /*=======================================
     *  generate U from the bulge
     *=======================================*/
    if (jobu == PlasmaVec){
        #if defined(ENABLE_TIMER)
        timeT  = PLASMA_Wtime();
        #endif
        /* compute T2 */
        plasma_static_call_8(plasma_pzlarft_blgtrd,
            int,                 MINMN,
            int,                 NB,
            int,                 Vblksiz,
            PLASMA_Complex64_t*, VQ2,
            PLASMA_Complex64_t*, TQ2,
            PLASMA_Complex64_t*, TAUQ2,
            PLASMA_sequence*,    sequence,
            PLASMA_request*,     request);
        #if defined(ENABLE_TIMER)
        plasma_dynamic_sync();
        timeT  = PLASMA_Wtime()-timeT;
        printf("  Finish compute TU2  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,                 MINMN,
            int,                 NB,
            int,                 MINMNE,
            int,                 Vblksiz,
            int,                 WANTZ,
            PLASMA_Complex64_t*, VQ2,
            PLASMA_Complex64_t*, TQ2,
            PLASMA_Complex64_t*, TAUQ2,
            PLASMA_Complex64_t*, U,
            int,                 LDU,
            PLASMA_sequence*,    sequence,
            PLASMA_request*,     request);
        #if defined(ENABLE_TIMER)
        plasma_dynamic_sync();
        timeaplQ2  = PLASMA_Wtime()-timeaplQ2;
        printf("  Finish compute U2  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( descU, U, NB, NB, LDU, ME, 0, 0, M, ME, sequence, request, plasma_desc_mat_free(&(descU)) );
            } else {
                plasma_ziplap2tile( descU, U, NB, NB, LDU, ME, 0, 0, M, ME, sequence, request);
            }
            #if defined(ENABLE_TIMER)
            timeconv1    = PLASMA_Wtime()-timeconv1;
            timeaplQ1   = PLASMA_Wtime();
            #endif
            /* Accumulate the transformations from the first stage */
            if(M<N){
                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(descU, descU.mb, 0, descU.m-descU.mb, descU.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_pzunmqr,
                    PLASMA_enum, PlasmaLeft,
                    PLASMA_enum, PlasmaNoTrans,
                    PLASMA_desc, descA,
                    PLASMA_desc, descU,
                    PLASMA_desc, descT,
                    PLASMA_sequence*, sequence,
                    PLASMA_request*, request);
            }
            #if defined(ENABLE_TIMER)
            plasma_dynamic_sync();
            timeaplQ1   = PLASMA_Wtime()-timeaplQ1;
            printf("  Finish compute U1  timing= %lf \n" ,timeaplQ1);
            timeconv2    = PLASMA_Wtime();
            #endif
            if ( PLASMA_TRANSLATION == PLASMA_OUTOFPLACE ) {
                plasma_zooptile2lap( descU, U, NB, NB, LDU, ME, sequence, request );
                plasma_dynamic_sync();
                plasma_desc_mat_free(&descU);
            } else {
                plasma_ziptile2lap( descU, U, NB, NB, LDU, ME, sequence, request );
                plasma_dynamic_sync();
            }
            #if defined(ENABLE_TIMER)
            plasma_dynamic_sync();
            timeconv2    = PLASMA_Wtime()-timeconv2;
            printf("  Finish convert  U  timing= %lf \n" ,timeconv1+timeconv2);
            #endif
        } /* END of ( NB < N ) */
    }
    /*=======================================
     *  END of calling computing U
     *=======================================*/
    /*=======================================
     *  generate VT from the bulge
     *=======================================*/
    if (jobvt == PlasmaVec){
        #if defined(ENABLE_TIMER)
        timeT  = PLASMA_Wtime();
        #endif
        /* compute T2 */
        plasma_static_call_8(plasma_pzlarft_blgtrd,
            int,                 MINMN,
            int,                 NB,
            int,                 Vblksiz,
            PLASMA_Complex64_t*, VP2,
            PLASMA_Complex64_t*, TP2,
            PLASMA_Complex64_t*, TAUP2,
            PLASMA_sequence*,    sequence,
            PLASMA_request*,     request);
        #if defined(ENABLE_TIMER)
        plasma_dynamic_sync();
        timeT  = PLASMA_Wtime()-timeT;
        printf("  Finish compute TV2  timing= %lf \n" ,timeT);
        timeaplQ2  = PLASMA_Wtime();
        #endif
        /* apply Q2 from Left */
        plasma_static_call_14(plasma_pzunmqr_blgtrd,
            PLASMA_enum,         PlasmaRight,
            PLASMA_enum,         PlasmaConjTrans,
            int,                 MINMN,
            int,                 NB,
            int,                 MINMNE,
            int,                 Vblksiz,
            int,                 WANTZ,
            PLASMA_Complex64_t*, VP2,
            PLASMA_Complex64_t*, TP2,
            PLASMA_Complex64_t*, TAUP2,
            PLASMA_Complex64_t*, VT,
            int,                 LDVT,
            PLASMA_sequence*,    sequence,
            PLASMA_request*,     request);
        #if defined(ENABLE_TIMER)
        plasma_dynamic_sync();
        timeaplQ2  = PLASMA_Wtime()-timeaplQ2;
        printf("  Finish compute V2 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( descVT, VT, NB, NB, LDVT, N, 0, 0, NE, N, sequence, request, plasma_desc_mat_free(&(descVT)) );
            } else {
                plasma_ziplap2tile( descVT, VT, NB, NB, LDVT, N, 0, 0, NE, N, sequence, request);
            }
            #if defined(ENABLE_TIMER)
            timeconv1    = PLASMA_Wtime()-timeconv1;
            timeaplQ1   = PLASMA_Wtime();
            #endif
            /* Accumulate the transformations from the first stage */
            if(M<N){
                plasma_dynamic_call_7(plasma_pzunmlq,
                    PLASMA_enum, PlasmaRight,
                    PLASMA_enum, PlasmaNoTrans,
                    PLASMA_desc, descA,
                    PLASMA_desc, descVT,
                    PLASMA_desc, descT,
                    PLASMA_sequence*, sequence,
                    PLASMA_request*, request);
            }
            else {
                plasma_dynamic_call_7(plasma_pzunmlq,
                    PLASMA_enum, PlasmaRight,
                    PLASMA_enum, PlasmaNoTrans,
                    PLASMA_desc, plasma_desc_submatrix(descA, 0, descA.nb, descA.m-descA.mb, descA.n-descA.nb),
                    PLASMA_desc, plasma_desc_submatrix(descVT,0, descVT.nb, descVT.m, descVT.n-descVT.nb),
                    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 V1  timing= %lf \n" ,timeaplQ1);
            timeconv2    = PLASMA_Wtime();
            #endif

            if ( PLASMA_TRANSLATION == PLASMA_OUTOFPLACE ) {
                plasma_zooptile2lap( descVT, VT, NB, NB, LDVT, N, sequence, request );
                plasma_dynamic_sync();
                plasma_desc_mat_free(&descVT);
            } else {
                plasma_ziptile2lap( descVT, VT, NB, NB, LDVT, N, sequence, request );
                plasma_dynamic_sync();
            }
            #if defined(ENABLE_TIMER)
            plasma_dynamic_sync();
            timeconv2    = PLASMA_Wtime()-timeconv2;
            printf("  Finish convert VT  timing= %lf \n" ,timeconv1+timeconv2);
            #endif
        } /* END of ( NB < N ) */
    }
    /*=======================================
     *  END of calling computing VT
     *=======================================*/
    #if defined(ENABLE_TIMER)
    timeall = PLASMA_Wtime()-timeall;
    printf("  Finish full eigenproblem threads %d  N %d  fraction %d  MINMNE %d  timeall= %lf \n", plasma->world_size, MINMN, fraction, MINMNE, timeall);
    #endif
    if( jobu  == PlasmaVec )
        plasma_shared_free(plasma, TQ2);
    if( jobvt == PlasmaVec )
        plasma_shared_free(plasma, TP2);
    plasma_shared_free(plasma, VQ2);
    plasma_shared_free(plasma, TAUQ2);
    plasma_shared_free(plasma, VP2);
    plasma_shared_free(plasma, TAUP2);
    plasma_shared_free(plasma, E);
    plasma_shared_free(plasma, AB);
    return PLASMA_SUCCESS;
}