text:
#include <stdio.h>
#include <stdarg.h>
#include <string.h>
#include <math.h>
#include <gbpLib.h>
#include <gbpRNG.h>
#include <gbpMCMC.h>
#include <gsl/gsl_linalg.h>
#include <gsl/gsl_fit.h>
#include <gsl/gsl_interp.h>

void read_MCMC_state(MCMC_info *MCMC) {
    char          filename_output_dir[SID_MAX_FILENAME_LENGTH];
    char          filename_chain_dir[SID_MAX_FILENAME_LENGTH];
    char          filename_results_dir[SID_MAX_FILENAME_LENGTH];
    char          filename_plots_dir[SID_MAX_FILENAME_LENGTH];
    char          filename_run[SID_MAX_FILENAME_LENGTH];
    char          filename_chain[SID_MAX_FILENAME_LENGTH];
    char          filename_chain_config[SID_MAX_FILENAME_LENGTH];
    char          filename_stats[SID_MAX_FILENAME_LENGTH];
    char          filename_coverage[SID_MAX_FILENAME_LENGTH];
    char          filename_chain_covariance[SID_MAX_FILENAME_LENGTH];
    char          filename_covariance[SID_MAX_FILENAME_LENGTH];
    char          filename_histograms[SID_MAX_FILENAME_LENGTH];
    char          filename_results[SID_MAX_FILENAME_LENGTH];
    char          filename_stop[SID_MAX_FILENAME_LENGTH];
    char          format_string[32];
    int           my_chain;
    int           i_P, i_DS, i_M, i_array;
    double *      V_read;
    FILE *        fp_run;
    FILE *        fp_chain;
    FILE *        fp_chain_config;
    FILE *        fp_stats;
    FILE *        fp_coverage;
    FILE *        fp_chain_covariance;
    FILE *        fp_covariance;
    FILE *        fp_histograms;
    FILE *        fp_results;
    FILE *        fp_stop;
    MCMC_DS_info *current_DS;

    set_MCMC_mode(MCMC, MCMC_MODE_DEFAULT);
    my_chain = MCMC->my_chain;

    SID_log("Reading MCMC state from {%s}...", SID_LOG_OPEN, MCMC->filename_output_dir);

    // Set directories
    sprintf(filename_output_dir, "%s/", MCMC->filename_output_dir);
    sprintf(filename_chain_dir, "%s/chains/", MCMC->filename_output_dir);
    sprintf(filename_results_dir, "%s/results/", MCMC->filename_output_dir);
    sprintf(filename_plots_dir, "%s/plots/", MCMC->filename_output_dir);
    // Set filenames
    sprintf(filename_run, "%s/run.dat", MCMC->filename_output_dir);
    sprintf(filename_chain, "%s/chain_trace_%06d.dat", filename_chain_dir, my_chain);
    sprintf(filename_chain_config, "%s/chain_config_%06d.dat", filename_chain_dir, my_chain);
    sprintf(filename_chain_covariance, "%s/chain_covariance_%06d.dat", filename_chain_dir, my_chain);
    sprintf(filename_stats, "%s/chain_stats_%06d.dat", filename_chain_dir, my_chain);
    sprintf(filename_coverage, "%s/coverage.dat", filename_results_dir);
    sprintf(filename_histograms, "%s/histograms.dat", filename_results_dir);
    sprintf(filename_covariance, "%s/covariance.dat", filename_results_dir);

    MCMC->map_P_to_M                 = NULL;
    MCMC->compute_MCMC_ln_likelihood = compute_MCMC_ln_likelihood_default;
    MCMC->params                     = NULL;
    MCMC->temperature                = 1.0;
    MCMC->n_P                        = 0;
    MCMC->n_thin                     = 1;
    MCMC->n_DS                       = 0;
    MCMC->n_M_total                  = 0;
    MCMC->n_arrays                   = 0;
    MCMC->n_M                        = NULL;
    MCMC->array                      = NULL;
    MCMC->V                          = NULL;
    MCMC->m                          = NULL;
    MCMC->b                          = NULL;
    MCMC->RNG                        = NULL;
    MCMC->flag_integrate_on          = GBP_TRUE;
    MCMC->flag_analysis_on           = GBP_TRUE;
    MCMC->first_map_call             = GBP_TRUE;
    MCMC->mode                       = MCMC_MODE_DEFAULT;
    MCMC->DS                         = NULL;
    MCMC->last                       = NULL;

    // Read/Write Header file
    if((fp_run = fopen(filename_run, "rb")) != NULL) {
        fp_run = fopen(filename_run, "rb");
        SID_fread_verify(MCMC->problem_name, sizeof(char), MCMC_NAME_SIZE, fp_run);
        SID_fread_verify(&(MCMC->n_chains), sizeof(int), 1, fp_run);
        SID_fread_verify(&(MCMC->n_avg), sizeof(int), 1, fp_run);
        SID_fread_verify(&(MCMC->flag_autocor_on), sizeof(int), 1, fp_run);
        SID_fread_verify(&(MCMC->flag_no_map_write), sizeof(int), 1, fp_run);
        SID_fread_verify(&(MCMC->n_P), sizeof(int), 1, fp_run);
        SID_log("Problem name    ={%s}", SID_LOG_COMMENT, MCMC->problem_name);
        SID_log("n_avg           ={%d}", SID_LOG_COMMENT, MCMC->n_avg);
        SID_log("flag_autocor_on ={%d}", SID_LOG_COMMENT, MCMC->flag_autocor_on);
        MCMC->P_names     = (char **)SID_malloc(sizeof(char *) * MCMC->n_P);
        MCMC->P_init      = (double *)SID_malloc(sizeof(double) * MCMC->n_P);
        MCMC->P_new       = (double *)SID_malloc(sizeof(double) * MCMC->n_P);
        MCMC->P_last      = (double *)SID_malloc(sizeof(double) * MCMC->n_P);
        MCMC->P_chain     = (double *)SID_malloc(sizeof(double) * MCMC->n_P);
        MCMC->P_limit_min = (double *)SID_malloc(sizeof(double) * MCMC->n_P);
        MCMC->P_limit_max = (double *)SID_malloc(sizeof(double) * MCMC->n_P);
        for(i_P = 0; i_P < MCMC->n_P; i_P++)
            MCMC->P_limit_min[i_P] = -DBL_MAX * 1e-3;
        for(i_P = 0; i_P < MCMC->n_P; i_P++)
            MCMC->P_limit_max[i_P] = DBL_MAX * 1e-3;

        SID_log("Parameters (name, initial_value,limit min,limit max):", SID_LOG_OPEN);
        MCMC->P_name_length = 0;
        for(i_P = 0; i_P < MCMC->n_P; i_P++) {
            MCMC->P_names[i_P] = (char *)SID_malloc(sizeof(char) * MCMC_NAME_SIZE);
            SID_fread_verify(MCMC->P_names[i_P], sizeof(char), MCMC_NAME_SIZE, fp_run);
            SID_fread_verify(&(MCMC->P_init[i_P]), sizeof(double), 1, fp_run);
            SID_fread_verify(&(MCMC->P_limit_min[i_P]), sizeof(double), 1, fp_run);
            SID_fread_verify(&(MCMC->P_limit_max[i_P]), sizeof(double), 1, fp_run);
            MCMC->P_name_length = GBP_MAX((size_t)(MCMC->P_name_length), strlen(MCMC->P_names[i_P]));
        }
        sprintf(MCMC->P_name_format, "%%-%ds", MCMC->P_name_length);
        sprintf(format_string, "%s %%13.6le %%13.6le %%13.6le", MCMC->P_name_format);
        for(i_P = 0; i_P < MCMC->n_P; i_P++)
            SID_log(format_string, SID_LOG_COMMENT, MCMC->P_names[i_P], MCMC->P_init[i_P], MCMC->P_limit_min[i_P], MCMC->P_limit_max[i_P]);
        SID_log(NULL, SID_LOG_CLOSE | SID_LOG_NOPRINT);
        SID_fread_verify(&(MCMC->n_arrays), sizeof(int), 1, fp_run);
        SID_log("n_arrays=%d", SID_LOG_OPEN, MCMC->n_arrays);
        MCMC->array      = (double **)SID_malloc(sizeof(double *) * MCMC->n_arrays);
        MCMC->array_name = (char **)SID_malloc(sizeof(char *) * MCMC->n_arrays);
        for(i_array = 0; i_array < MCMC->n_arrays; i_array++) {
            MCMC->array[i_array]      = (double *)SID_malloc(sizeof(double) * MCMC->n_P);
            MCMC->array_name[i_array] = (char *)SID_malloc(sizeof(char) * MCMC_NAME_SIZE);
            SID_fread_verify(MCMC->array_name[i_array], sizeof(char), MCMC_NAME_SIZE, fp_run);
            SID_fread_verify(MCMC->array[i_array], sizeof(double), MCMC->n_P, fp_run);
            SID_log("array #%03d name ={%s}", SID_LOG_COMMENT, i_array, MCMC->array_name[i_array]);
        }
        SID_log(NULL, SID_LOG_CLOSE | SID_LOG_NOPRINT);
        SID_fread_verify(&(MCMC->n_DS), sizeof(int), 1, fp_run);
        SID_log("Reading %d datasets...", SID_LOG_OPEN, MCMC->n_DS);
        for(i_DS = 0; i_DS < MCMC->n_DS; i_DS++) {
            SID_log("Dataset #%03d:", SID_LOG_OPEN, i_DS);
            current_DS = (MCMC_DS_info *)SID_malloc(sizeof(MCMC_DS_info));
            SID_fread_verify(current_DS->name, sizeof(char), MCMC_NAME_SIZE, fp_run);
            SID_fread_verify(&(current_DS->n_M), sizeof(int), 1, fp_run);
            MCMC->n_M_total += current_DS->n_M;
            SID_log("name    ={%s}", SID_LOG_COMMENT, current_DS->name);
            SID_log("n_M     =%d", SID_LOG_COMMENT, current_DS->n_M);
            current_DS->M_target  = (double *)SID_malloc(sizeof(double) * current_DS->n_M);
            current_DS->dM_target = (double *)SID_malloc(sizeof(double) * current_DS->n_M);
            current_DS->params    = NULL;
            SID_fread_verify(current_DS->M_target, sizeof(double), current_DS->n_M, fp_run);
            SID_fread_verify(current_DS->dM_target, sizeof(double), current_DS->n_M, fp_run);
            SID_fread_verify(&(current_DS->n_arrays), sizeof(int), 1, fp_run);
            SID_log("n_arrays=%d", SID_LOG_OPEN, current_DS->n_arrays);
            current_DS->array      = (double **)SID_malloc(sizeof(double *) * current_DS->n_arrays);
            current_DS->array_name = (char **)SID_malloc(sizeof(char *) * current_DS->n_arrays);
            for(i_array = 0; i_array < current_DS->n_arrays; i_array++) {
                current_DS->array[i_array]      = (double *)SID_malloc(sizeof(double) * current_DS->n_M);
                current_DS->array_name[i_array] = (char *)SID_malloc(sizeof(char) * MCMC_NAME_SIZE);
                SID_fread_verify(current_DS->array_name[i_array], sizeof(char), MCMC_NAME_SIZE, fp_run);
                SID_fread_verify(current_DS->array[i_array], sizeof(double), current_DS->n_M, fp_run);
                SID_log("array #%03d name={%s}", SID_LOG_COMMENT, i_array, current_DS->array_name[i_array]);
            }
            SID_log(NULL, SID_LOG_CLOSE | SID_LOG_NOPRINT);
            current_DS->next = NULL;
            if(MCMC->DS == NULL)
                MCMC->DS = current_DS;
            else
                MCMC->last->next = current_DS;
            MCMC->last = current_DS;
            SID_log(NULL, SID_LOG_CLOSE | SID_LOG_NOPRINT);
        }
        SID_log("Done.", SID_LOG_CLOSE);
        fclose(fp_run);
    }

    // ... fetch the number of intervals that have already been computed ...
    fp_chain_config = fopen(filename_chain_config, "rb");
    V_read          = (double *)SID_malloc(sizeof(double) * MCMC->n_P * MCMC->n_P);
    SID_fread_verify(&(MCMC->n_iterations), sizeof(int), 1, fp_chain_config);
    SID_fread_verify(&(MCMC->n_iterations_burn), sizeof(int), 1, fp_chain_config);

    // ... fetch the temperature and covariance matrix that was being used
    SID_fread_verify(&(MCMC->temperature), sizeof(double), 1, fp_chain_config);
    SID_fread_verify(V_read, sizeof(double), MCMC->n_P * MCMC->n_P, fp_chain_config);
    set_MCMC_covariance(MCMC, V_read);
    SID_free(SID_FARG V_read);

    // Initialize dataset arrays
    init_MCMC_arrays(MCMC);

    SID_log("# burn  iterations = %d", SID_LOG_COMMENT, MCMC->n_iterations_burn);
    SID_log("# total iterations = %d", SID_LOG_COMMENT, MCMC->n_iterations);
    SID_log("Temperature        = %le", SID_LOG_COMMENT, MCMC->temperature);
    fclose(fp_chain_config);

    SID_log("Done.", SID_LOG_CLOSE);
}

ext: c

max_stars_repo_name: gbpoole/gbpCode

max_issues_repo_name: gbpoole/gbpCode

max_forks_repo_name: gbpoole/gbpCode


text:
/**
 *
 * @file core_dgbtype2cb.c
 *
 *  PLASMA core_blas kernel
 *  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 2012-12-15
 * @generated d Tue Jan  7 11:44:50 2014
 *
 **/
#include <lapacke.h>
#include "common.h"

//#define AU(m,n) &(A[(m) + LDA*(n)])
//#define AL(m,n) &(A[(m) + LDA*(n)])
#define AL(m_, n_) (A + NB + LDA * (n_) + ((m_)-(n_)))
#define AU(m_, n_) (A + NB + LDA * (n_) + ((m_)-(n_)+NB))
#define VQ(m)     (VQ + (m))
#define VP(m)     (VP + (m))
#define TAUQ(m)   (TAUQ + (m))
#define TAUP(m)   (TAUP + (m))

/***************************************************************************//**
 *
 * @ingroup CORE_double
 *
 *  CORE_dgbtype2cb is a kernel that will operate on a region (triangle) of data
 *  bounded by st and ed. This kernel apply the right update remaining from the
 *  type1 and this later will create a bulge so it eliminate the first column of
 *  the created bulge and do the corresponding Left update.
 *
 *  All detail are available on technical report or SC11 paper.
 *  Azzam Haidar, Hatem Ltaief, and Jack Dongarra. 2011.
 *  Parallel reduction to condensed forms for symmetric eigenvalue problems
 *  using aggregated fine-grained and memory-aware kernels. In Proceedings
 *  of 2011 International Conference for High Performance Computing,
 *  Networking, Storage and Analysis (SC '11). ACM, New York, NY, USA, ,
 *  Article 8 , 11 pages.
 *  http://doi.acm.org/10.1145/2063384.2063394
 *
 *******************************************************************************
 *
 * @param[in] N
 *          The order of the matrix A.
 *
 * @param[in] NB
 *          The size of the band.
 *
 * @param[in, out] A
 *          A pointer to the matrix A of size (3*NB+1)-by-N.
 *
 * @param[in] LDA
 *          The leading dimension of the matrix A. LDA >= max(1,3*NB+1)
 *
 * @param[in, out] V
 *          double array, dimension N if eigenvalue only
 *          requested or (LDV*blkcnt*Vblksiz) if Eigenvectors requested
 *          The Householder reflectors of the previous type 1 are used here
 *          to continue update then new one are generated to eliminate the
 *          bulge and stored in this array.
 *
 * @param[in, out] TAU
 *          double array, dimension (N).
 *          The scalar factors of the Householder reflectors of the previous
 *          type 1 are used here to continue update then new one are generated
 *          to eliminate the bulge and stored in this array.
 *
 * @param[in] st
 *          A pointer to the start index where this kernel will operate.
 *
 * @param[in] ed
 *          A pointer to the end index where this kernel will operate.
 *
 * @param[in] sweep
 *          The sweep number that is eliminated. it serve to calculate the
 *          pointer to the position where to store the Vs and Ts.
 *
 * @param[in] Vblksiz
 *          constant which correspond to the blocking used when applying the Vs.
 *          it serve to calculate the pointer to the position where to store the
 *          Vs and Ts.
 *
 * @param[in] WANTZ
 *          constant which indicate if Eigenvalue are requested or both
 *          Eigenvalue/Eigenvectors.
 *
 * @param[in] WORK
 *          Workspace of size nb.
 *
 *******************************************************************************
 *
 * @return
 *          \retval PLASMA_SUCCESS successful exit
 *          \retval <0 if -i, the i-th argument had an illegal value
 *
 ******************************************************************************/
/***************************************************************************
 *          TYPE 2-BAND-bidiag Lower/Upper columnwise-Householder
 ***************************************************************************/
void
CORE_dgbtype2cb(PLASMA_enum uplo, int N, int NB,
                double *A, int LDA,
                double *VQ, double *TAUQ,
                double *VP, double *TAUP,
                int st, int ed, int sweep, int Vblksiz, int WANTZ,
                double *WORK)
{
    double ctmp;
    int i, J1, J2, len, lem, LDX;
    int blkid, vpos, taupos, tpos;

    LDX = LDA-1;
    J1  = ed+1;
    J2  = min(ed+NB,N-1);
    lem = ed-st+1;
    len = J2-J1+1;

    if( uplo == PlasmaUpper ) {
        /* ========================
        *       UPPER CASE
        * ========================*/
        if( len > 0 ) {
            if( WANTZ == 0 ) {
                vpos   = ((sweep+1)%2)*N + st;
                taupos = ((sweep+1)%2)*N + st;
            } else {
                findVTpos(N, NB, Vblksiz, sweep, st,
                          &vpos, &taupos, &tpos, &blkid);
            }
            /* Apply remaining Left commming from type1/3_upper */
            ctmp = (*TAUQ(taupos));
            LAPACKE_dlarfx_work(LAPACK_COL_MAJOR, lapack_const(PlasmaLeft),
                                lem, len, VQ(vpos), ctmp, AU(st, J1), LDX, WORK);
        }

        if( len > 1 ) {
            if( WANTZ == 0 ) {
                vpos   = ((sweep+1)%2)*N + J1;
                taupos = ((sweep+1)%2)*N + J1;
            } else {
                findVTpos(N,NB,Vblksiz,sweep,J1, &vpos, &taupos, &tpos, &blkid);
            }

            /* Remove the top row of the created bulge */
            *VP(vpos) = 1.;
            for(i=1; i<len; i++){
                *VP(vpos+i)     = (*AU(st, J1+i));
                *AU(st, J1+i) = 0.;
            }
            /* Eliminate the row  at st */
            ctmp = (*AU(st, J1));
            LAPACKE_dlarfg_work(len, &ctmp, VP(vpos+1), 1, TAUP(taupos) );
            *AU(st, J1) = ctmp;
            /*
             * Apply Right on A(J1:J2,st+1:ed)
             * We decrease len because we start at row st+1 instead of st.
             * row st is the row that has been revomved;
             */
            lem = lem-1;
            ctmp = *TAUP(taupos);
            LAPACKE_dlarfx_work(LAPACK_COL_MAJOR, lapack_const(PlasmaRight),
                                lem, len, VP(vpos), ctmp, AU(st+1, J1), LDX, WORK);
        }
    }else{
        /* ========================
         *       LOWER CASE
         * ========================*/
        if( len > 0 ) {
            if( WANTZ == 0 ) {
                vpos   = ((sweep+1)%2)*N + st;
                taupos = ((sweep+1)%2)*N + st;
            } else {
                findVTpos(N, NB, Vblksiz, sweep, st,
                          &vpos, &taupos, &tpos, &blkid);
            }
            /* Apply remaining Right commming from type1/3_lower */
            ctmp = (*TAUP(taupos));
            LAPACKE_dlarfx_work(LAPACK_COL_MAJOR, lapack_const(PlasmaRight),
                                len, lem, VP(vpos), ctmp, AL(J1, st), LDX, WORK);
        }
        if( len > 1 ) {
            if( WANTZ == 0 ) {
                vpos   = ((sweep+1)%2)*N + J1;
                taupos = ((sweep+1)%2)*N + J1;
            } else {
                findVTpos(N,NB,Vblksiz,sweep,J1, &vpos, &taupos, &tpos, &blkid);
            }

            /* Remove the first column of the created bulge */
            *VQ(vpos) = 1.;
            memcpy(VQ(vpos+1), AL(J1+1, st), (len-1)*sizeof(double));
            memset(AL(J1+1, st), 0, (len-1)*sizeof(double));
            /* Eliminate the col  at st */
            LAPACKE_dlarfg_work(len, AL(J1, st), VQ(vpos+1), 1, TAUQ(taupos) );
            /*
             * Apply left on A(J1:J2,st+1:ed)
             * We decrease len because we start at col st+1 instead of st.
             * col st is the col that has been revomved;
             */
            lem = lem-1;
            ctmp = (*TAUQ(taupos));
            LAPACKE_dlarfx_work(LAPACK_COL_MAJOR, lapack_const(PlasmaLeft),
                                len, lem, VQ(vpos), ctmp, AL(J1, st+1), LDX, WORK);
        }
    }
    /* end of uplo case */
    return;
}
/***************************************************************************/
#undef AU
#undef AL
#undef VQ
#undef VP
#undef TAUQ
#undef TAUP


text:
/******************************** Inclusions. ********************************/
#include <stdio.h>
#include <stdlib.h>
#include <math.h>
#include <time.h>
#include <gsl/gsl_rng.h>
#include <gsl/gsl_randist.h>

/******************************** Definitions. *******************************/
#define MAX(x, y) ((x)>(y) ? (x) : (y))
#define MIN(x, y) ((x)<(y) ? (x) : (y))

ext: h

max_stars_repo_name: mgastner/impact-of-hypocrisy

max_issues_repo_name: mgastner/impact-of-hypocrisy

max_forks_repo_name: mgastner/impact-of-hypocrisy


text:
// Copyright (c) 2021 Stig Rune Sellevag
//
// This file is distributed under the MIT License. See the accompanying file
// LICENSE.txt or http://www.opensource.org/licenses/mit-license.php for terms
// and conditions.

#ifndef SCILIB_LINALG_LSTSQ_H
#define SCILIB_LINALG_LSTSQ_H

#ifdef USE_MKL
#include <mkl.h>
#else
#include <lapacke.h>
#endif

#include <scilib/mdarray.h>
#include <scilib/linalg_impl/lapack_types.h>
#include <experimental/mdspan>
#include <exception>
#include <algorithm>
#include <type_traits>

namespace Sci {
namespace Linalg {

// Compute the minimum norm-solution to a real linear least squares problem.
template <class Layout>
inline void lstsq(Sci::Matrix_view<double, Layout> a,
                  Sci::Matrix_view<double, Layout> b)
{
    namespace stdex = std::experimental;

    static_assert(a.is_contiguous());
    static_assert(b.is_contiguous());

    BLAS_INT m = static_cast<BLAS_INT>(a.extent(0));
    BLAS_INT n = static_cast<BLAS_INT>(a.extent(1));
    BLAS_INT nrhs = static_cast<BLAS_INT>(b.extent(1));
    BLAS_INT rank;

    double rcond = -1.0;                           // use machine epsilon
    Sci::Vector<double, Layout> s(std::min(m, n)); // singular values of a

    auto matrix_layout = LAPACK_ROW_MAJOR;
    BLAS_INT lda = n;
    BLAS_INT ldb = nrhs;

    if constexpr (std::is_same_v<Layout, stdex::layout_left>) {
        matrix_layout = LAPACK_COL_MAJOR;
        lda = m;
        ldb = n;
    }
    BLAS_INT info = LAPACKE_dgelsd(matrix_layout, m, n, nrhs, a.data(), lda,
                                   b.data(), ldb, s.data(), rcond, &rank);
    if (info != 0) {
        throw std::runtime_error("dgelsd failed");
    }
}

template <class Layout, class Allocator>
inline void lstsq(Sci::Matrix<double, Layout, Allocator>& a,
                  Sci::Matrix<double, Layout, Allocator>& b)
{
    lstsq(a.view(), b.view());
}

} // namespace Linalg
} // namespace Sci

#endif // SCILIB_LINALG_LSTSQ_H

max_stars_repo_name: stigrs/scilib

max_issues_repo_name: stigrs/scilib

max_forks_repo_name: stigrs/scilib


text:
/* eigen/nonsymmv.c
 *
 * Copyright (C) 2006 Patrick Alken
 *
 * This program is free software; you can redistribute it and/or modify
 * it under the terms of the GNU General Public License as published by
 * the Free Software Foundation; either version 2 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 <stdlib.h>
#include <math.h>

#include <gsl/gsl_complex.h>
#include <gsl/gsl_complex_math.h>
#include <gsl/gsl_eigen.h>
#include <gsl/gsl_linalg.h>
#include <gsl/gsl_math.h>
#include <gsl/gsl_blas.h>
#include <gsl/gsl_vector.h>
#include <gsl/gsl_vector_complex.h>
#include <gsl/gsl_matrix.h>

/*
 * This module computes the eigenvalues and eigenvectors of a real
 * nonsymmetric matrix.
 *
 * This file contains routines based on original code from LAPACK
 * which is distributed under the modified BSD license. The LAPACK
 * routines used are DTREVC and DLALN2.
 */

#define GSL_NONSYMMV_SMLNUM (2.0 * GSL_DBL_MIN)
#define GSL_NONSYMMV_BIGNUM ((1.0 - GSL_DBL_EPSILON) / GSL_NONSYMMV_SMLNUM)

static void nonsymmv_get_right_eigenvectors(gsl_matrix *T, gsl_matrix *Z,
                                            gsl_vector_complex *eval,
                                            gsl_matrix_complex *evec,
                                            gsl_eigen_nonsymmv_workspace *w);
static inline void nonsymmv_solve_equation(gsl_matrix *A, double z,
                                           gsl_vector *b, gsl_vector *x,
                                           double *s, double *xnorm,
                                           double smin);
static inline void nonsymmv_solve_equation_z(gsl_matrix *A, gsl_complex *z,
                                             gsl_vector_complex *b,
                                             gsl_vector_complex *x,
                                             double *s, double *xnorm,
                                             double smin);
static void nonsymmv_normalize_eigenvectors(gsl_vector_complex *eval,
                                            gsl_matrix_complex *evec);

/*
gsl_eigen_nonsymmv_alloc()

Allocate a workspace for solving the nonsymmetric eigenvalue problem.
The size of this workspace is O(5n).

Inputs: n - size of matrices

Return: pointer to workspace
*/

gsl_eigen_nonsymmv_workspace *
gsl_eigen_nonsymmv_alloc(const size_t n)
{
  gsl_eigen_nonsymmv_workspace *w;

  if (n == 0)
    {
      GSL_ERROR_NULL ("matrix dimension must be positive integer",
                      GSL_EINVAL);
    }

  w = (gsl_eigen_nonsymmv_workspace *)
      malloc (sizeof (gsl_eigen_nonsymmv_workspace));

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

  w->size = n;
  w->Z = NULL;
  w->nonsymm_workspace_p = gsl_eigen_nonsymm_alloc(n);

  if (w->nonsymm_workspace_p == 0)
    {
      GSL_ERROR_NULL ("failed to allocate space for nonsymm workspace", GSL_ENOMEM);
    }

  /*
   * set parameters to compute the full Schur form T and balance
   * the matrices
   */
  gsl_eigen_nonsymm_params(1, 1, w->nonsymm_workspace_p);

  w->work = gsl_vector_alloc(n);
  w->work2 = gsl_vector_alloc(n);
  w->work3 = gsl_vector_alloc(n);
  if (w->work == 0 || w->work2 == 0 || w->work3 == 0)
    {
      GSL_ERROR_NULL ("failed to allocate space for nonsymmv additional workspace", GSL_ENOMEM);
    }

  return (w);
} /* gsl_eigen_nonsymmv_alloc() */

/*
gsl_eigen_nonsymmv_free()
  Free workspace w
*/

void
gsl_eigen_nonsymmv_free (gsl_eigen_nonsymmv_workspace * w)
{
  gsl_eigen_nonsymm_free(w->nonsymm_workspace_p);
  gsl_vector_free(w->work);
  gsl_vector_free(w->work2);
  gsl_vector_free(w->work3);

  free(w);
} /* gsl_eigen_nonsymmv_free() */

/*
gsl_eigen_nonsymmv()

Solve the nonsymmetric eigensystem problem

A x = \lambda x

for the eigenvalues \lambda and right eigenvectors x

Inputs: A    - general real matrix
        eval - where to store eigenvalues
        evec - where to store eigenvectors
        w    - workspace

Return: success or error
*/

int
gsl_eigen_nonsymmv (gsl_matrix * A, gsl_vector_complex * eval,
                    gsl_matrix_complex * evec,
                    gsl_eigen_nonsymmv_workspace * w)
{
  const size_t N = A->size1;

  /* check matrix and vector sizes */

  if (N != A->size2)
    {
      GSL_ERROR ("matrix must be square to compute eigenvalues", GSL_ENOTSQR);
    }
  else if (eval->size != N)
    {
      GSL_ERROR ("eigenvalue vector must match matrix size", GSL_EBADLEN);
    }
  else if (evec->size1 != evec->size2)
    {
      GSL_ERROR ("eigenvector matrix must be square", GSL_ENOTSQR);
    }
  else if (evec->size1 != N)
    {
      GSL_ERROR ("eigenvector matrix has wrong size", GSL_EBADLEN);
    }
  else
    {
      int s;
      gsl_matrix Z;

      /*
       * We need a place to store the Schur vectors, so we will
       * treat evec as a real matrix and store them in the left
       * half - the factor of 2 in the tda corresponds to the
       * complex multiplicity
       */
      Z.size1 = N;
      Z.size2 = N;
      Z.tda = 2 * N;
      Z.data = evec->data;
      Z.block = 0;
      Z.owner = 0;

      /* compute eigenvalues, Schur form, and Schur vectors */
      s = gsl_eigen_nonsymm_Z(A, eval, &Z, w->nonsymm_workspace_p);

      if (w->Z)
        {
          /*
           * save the Schur vectors in user supplied matrix, since
           * they will be destroyed when computing eigenvectors
           */
          gsl_matrix_memcpy(w->Z, &Z);
        }

      /* only compute eigenvectors if we found all eigenvalues */
      if (s == GSL_SUCCESS)
        {
          /* compute eigenvectors */
          nonsymmv_get_right_eigenvectors(A, &Z, eval, evec, w);

          /* normalize so that Euclidean norm is 1 */
          nonsymmv_normalize_eigenvectors(eval, evec);
        }

      return s;
    }
} /* gsl_eigen_nonsymmv() */

/*
gsl_eigen_nonsymmv_Z()
  Compute eigenvalues and eigenvectors of a real nonsymmetric matrix
and also save the Schur vectors. See comments in gsl_eigen_nonsymm_Z
for more information.

Inputs: A    - real nonsymmetric matrix
        eval - where to store eigenvalues
        evec - where to store eigenvectors
        Z    - where to store Schur vectors
        w    - nonsymmv workspace

Return: success or error
*/

int
gsl_eigen_nonsymmv_Z (gsl_matrix * A, gsl_vector_complex * eval,
                      gsl_matrix_complex * evec, gsl_matrix * Z,
                      gsl_eigen_nonsymmv_workspace * w)
{
  /* check matrix and vector sizes */

  if (A->size1 != A->size2)
    {
      GSL_ERROR ("matrix must be square to compute eigenvalues/eigenvectors", GSL_ENOTSQR);
    }
  else if (eval->size != A->size1)
    {
      GSL_ERROR ("eigenvalue vector must match matrix size", GSL_EBADLEN);
    }
  else if (evec->size1 != evec->size2)
    {
      GSL_ERROR ("eigenvector matrix must be square", GSL_ENOTSQR);
    }
  else if (evec->size1 != A->size1)
    {
      GSL_ERROR ("eigenvector matrix has wrong size", GSL_EBADLEN);
    }
  else if ((Z->size1 != Z->size2) || (Z->size1 != A->size1))
    {
      GSL_ERROR ("Z matrix has wrong dimensions", GSL_EBADLEN);
    }
  else
    {
      int s;

      w->Z = Z;

      s = gsl_eigen_nonsymmv(A, eval, evec, w);

      w->Z = NULL;

      return s;
    }
} /* gsl_eigen_nonsymmv_Z() */

/********************************************
 *           INTERNAL ROUTINES              *
 ********************************************/

/*
nonsymmv_get_right_eigenvectors()
  Compute the right eigenvectors of the Schur form T and then
backtransform them using the Schur vectors to get right eigenvectors of
the original matrix.

Inputs: T    - Schur form
        Z    - Schur vectors
        eval - where to store eigenvalues (to ensure that the
               correct eigenvalue is stored in the same position
               as the eigenvectors)
        evec - where to store eigenvectors
        w    - nonsymmv workspace

Return: none

Notes: 1) based on LAPACK routine DTREVC - the algorithm used is
          backsubstitution on the upper quasi triangular system T
          followed by backtransformation by Z to get vectors of the
          original matrix.

       2) The Schur vectors in Z are destroyed and replaced with
          eigenvectors stored with the same storage scheme as DTREVC.
          The eigenvectors are also stored in 'evec'

       3) The matrix T is unchanged on output

       4) Each eigenvector is normalized so that the element of
          largest magnitude has magnitude 1; here the magnitude of
          a complex number (x,y) is taken to be |x| + |y|
*/

static void
nonsymmv_get_right_eigenvectors(gsl_matrix *T, gsl_matrix *Z,
                                gsl_vector_complex *eval,
                                gsl_matrix_complex *evec,
                                gsl_eigen_nonsymmv_workspace *w)
{
  const size_t N = T->size1;
  const double smlnum = GSL_DBL_MIN * N / GSL_DBL_EPSILON;
  const double bignum = (1.0 - GSL_DBL_EPSILON) / smlnum;
  int i;              /* looping */
  size_t iu,          /* looping */
         ju,
         ii;
  gsl_complex lambda; /* current eigenvalue */
  double lambda_re,   /* Re(lambda) */
         lambda_im;   /* Im(lambda) */
  gsl_matrix_view Tv, /* temporary views */
                  Zv;
  gsl_vector_view y,  /* temporary views */
                  y2,
                  ev,
                  ev2;
  double dat[4],      /* scratch arrays */
         dat_X[4];
  double scale;       /* scale factor */
  double xnorm;       /* |X| */
  gsl_vector_complex_view ecol, /* column of evec */
                          ecol2;
  int complex_pair;   /* complex eigenvalue pair? */
  double smin;

  /*
   * Compute 1-norm of each column of upper triangular part of T
   * to control overflow in triangular solver
   */

  gsl_vector_set(w->work3, 0, 0.0);
  for (ju = 1; ju < N; ++ju)
    {
      gsl_vector_set(w->work3, ju, 0.0);
      for (iu = 0; iu < ju; ++iu)
        {
          gsl_vector_set(w->work3, ju,
                         gsl_vector_get(w->work3, ju) +
                         fabs(gsl_matrix_get(T, iu, ju)));
        }
    }

  for (i = (int) N - 1; i >= 0; --i)
    {
      iu = (size_t) i;

      /* get current eigenvalue and store it in lambda */
      lambda_re = gsl_matrix_get(T, iu, iu);

      if (iu != 0 && gsl_matrix_get(T, iu, iu - 1) != 0.0)
        {
          lambda_im = sqrt(fabs(gsl_matrix_get(T, iu, iu - 1))) *
                      sqrt(fabs(gsl_matrix_get(T, iu - 1, iu)));
        }
      else
        {
          lambda_im = 0.0;
        }

      GSL_SET_COMPLEX(&lambda, lambda_re, lambda_im);

      smin = GSL_MAX(GSL_DBL_EPSILON * (fabs(lambda_re) + fabs(lambda_im)),
                     smlnum);
      smin = GSL_MAX(smin, GSL_NONSYMMV_SMLNUM);

      if (lambda_im == 0.0)
        {
          int k, l;
          gsl_vector_view bv, xv;

          /* real eigenvector */

          /*
           * The ordering of eigenvalues in 'eval' is arbitrary and
           * does not necessarily follow the Schur form T, so store
           * lambda in the right slot in eval to ensure it corresponds
           * to the eigenvector we are about to compute
           */
          gsl_vector_complex_set(eval, iu, lambda);

          /*
           * We need to solve the system:
           *
           * (T(1:iu-1, 1:iu-1) - lambda*I)*X = -T(1:iu-1,iu)
           */

          /* construct right hand side */
          for (k = 0; k < i; ++k)
            {
              gsl_vector_set(w->work,
                             (size_t) k,
                             -gsl_matrix_get(T, (size_t) k, iu));
            }

          gsl_vector_set(w->work, iu, 1.0);

          for (l = i - 1; l >= 0; --l)
            {
              size_t lu = (size_t) l;

              if (lu == 0)
                complex_pair = 0;
              else
                complex_pair = gsl_matrix_get(T, lu, lu - 1) != 0.0;

              if (!complex_pair)
                {
                  double x;

                  /*
                   * 1-by-1 diagonal block - solve the system:
                   *
                   * (T_{ll} - lambda)*x = -T_{l(iu)}
                   */

                  Tv = gsl_matrix_submatrix(T, lu, lu, 1, 1);
                  bv = gsl_vector_view_array(dat, 1);
                  gsl_vector_set(&bv.vector, 0,
                                 gsl_vector_get(w->work, lu));
                  xv = gsl_vector_view_array(dat_X, 1);

                  nonsymmv_solve_equation(&Tv.matrix,
                                          lambda_re,
                                          &bv.vector,
                                          &xv.vector,
                                          &scale,
                                          &xnorm,
                                          smin);

                  /* scale x to avoid overflow */
                  x = gsl_vector_get(&xv.vector, 0);
                  if (xnorm > 1.0)
                    {
                      if (gsl_vector_get(w->work3, lu) > bignum / xnorm)
                        {
                          x /= xnorm;
                          scale /= xnorm;
                        }
                    }

                  if (scale != 1.0)
                    {
                      gsl_vector_view wv;

                      wv = gsl_vector_subvector(w->work, 0, iu + 1);
                      gsl_blas_dscal(scale, &wv.vector);
                    }

                  gsl_vector_set(w->work, lu, x);

                  if (lu > 0)
                    {
                      gsl_vector_view v1, v2;

                      /* update right hand side */

                      v1 = gsl_matrix_column(T, lu);
                      v1 = gsl_vector_subvector(&v1.vector, 0, lu);

                      v2 = gsl_vector_subvector(w->work, 0, lu);

                      gsl_blas_daxpy(-x, &v1.vector, &v2.vector);
                    } /* if (l > 0) */
                } /* if (!complex_pair) */
              else
                {
                  double x11, x21;

                  /*
                   * 2-by-2 diagonal block
                   */

                  Tv = gsl_matrix_submatrix(T, lu - 1, lu - 1, 2, 2);
                  bv = gsl_vector_view_array(dat, 2);
                  gsl_vector_set(&bv.vector, 0,
                                 gsl_vector_get(w->work, lu - 1));
                  gsl_vector_set(&bv.vector, 1,
                                 gsl_vector_get(w->work, lu));
                  xv = gsl_vector_view_array(dat_X, 2);

                  nonsymmv_solve_equation(&Tv.matrix,
                                          lambda_re,
                                          &bv.vector,
                                          &xv.vector,
                                          &scale,
                                          &xnorm,
                                          smin);

                  /* scale X(1,1) and X(2,1) to avoid overflow */
                  x11 = gsl_vector_get(&xv.vector, 0);
                  x21 = gsl_vector_get(&xv.vector, 1);

                  if (xnorm > 1.0)
                    {
                      double beta;

                      beta = GSL_MAX(gsl_vector_get(w->work3, lu - 1),
                                     gsl_vector_get(w->work3, lu));
                      if (beta > bignum / xnorm)
                        {
                          x11 /= xnorm;
                          x21 /= xnorm;
                          scale /= xnorm;
                        }
                    }

                  /* scale if necessary */
                  if (scale != 1.0)
                    {
                      gsl_vector_view wv;

                      wv = gsl_vector_subvector(w->work, 0, iu + 1);
                      gsl_blas_dscal(scale, &wv.vector);
                    }

                  gsl_vector_set(w->work, lu - 1, x11);
                  gsl_vector_set(w->work, lu, x21);

                  /* update right hand side */
                  if (lu > 1)
                    {
                      gsl_vector_view v1, v2;

                      v1 = gsl_matrix_column(T, lu - 1);
                      v1 = gsl_vector_subvector(&v1.vector, 0, lu - 1);
                      v2 = gsl_vector_subvector(w->work, 0, lu - 1);
                      gsl_blas_daxpy(-x11, &v1.vector, &v2.vector);

                      v1 = gsl_matrix_column(T, lu);
                      v1 = gsl_vector_subvector(&v1.vector, 0, lu - 1);
                      gsl_blas_daxpy(-x21, &v1.vector, &v2.vector);
                    }

                  --l;
                } /* if (complex_pair) */
            } /* for (l = i - 1; l >= 0; --l) */

          /*
           * At this point, w->work is an eigenvector of the
           * Schur form T. To get an eigenvector of the original
           * matrix, we multiply on the left by Z, the matrix of
           * Schur vectors
           */

          ecol = gsl_matrix_complex_column(evec, iu);
          y = gsl_matrix_column(Z, iu);

          if (iu > 0)
            {
              gsl_vector_view x;

              Zv = gsl_matrix_submatrix(Z, 0, 0, N, iu);

              x = gsl_vector_subvector(w->work, 0, iu);

              /* compute Z * w->work and store it in Z(:,iu) */
              gsl_blas_dgemv(CblasNoTrans,
                             1.0,
                             &Zv.matrix,
                             &x.vector,
                             gsl_vector_get(w->work, iu),
                             &y.vector);
            } /* if (iu > 0) */

          /* store eigenvector into evec */

          ev = gsl_vector_complex_real(&ecol.vector);
          ev2 = gsl_vector_complex_imag(&ecol.vector);

          scale = 0.0;
          for (ii = 0; ii < N; ++ii)
            {
              double a = gsl_vector_get(&y.vector, ii);

              /* store real part of eigenvector */
              gsl_vector_set(&ev.vector, ii, a);

              /* set imaginary part to 0 */
              gsl_vector_set(&ev2.vector, ii, 0.0);

              if (fabs(a) > scale)
                scale = fabs(a);
            }

          if (scale != 0.0)
            scale = 1.0 / scale;

          /* scale by magnitude of largest element */
          gsl_blas_dscal(scale, &ev.vector);
        } /* if (GSL_IMAG(lambda) == 0.0) */
      else
        {
          gsl_vector_complex_view bv, xv;
          size_t k;
          int l;
          gsl_complex lambda2;

          /* complex eigenvector */

          /*
           * Store the complex conjugate eigenvalues in the right
           * slots in eval
           */
          GSL_SET_REAL(&lambda2, GSL_REAL(lambda));
          GSL_SET_IMAG(&lambda2, -GSL_IMAG(lambda));
          gsl_vector_complex_set(eval, iu - 1, lambda);
          gsl_vector_complex_set(eval, iu, lambda2);

          /*
           * First solve:
           *
           * [ T(i:i+1,i:i+1) - lambda*I ] * X = 0
           */

          if (fabs(gsl_matrix_get(T, iu - 1, iu)) >=
              fabs(gsl_matrix_get(T, iu, iu - 1)))
            {
              gsl_vector_set(w->work, iu - 1, 1.0);
              gsl_vector_set(w->work2, iu,
                             lambda_im / gsl_matrix_get(T, iu - 1, iu));
            }
          else
            {
              gsl_vector_set(w->work, iu - 1,
                             -lambda_im / gsl_matrix_get(T, iu, iu - 1));
              gsl_vector_set(w->work2, iu, 1.0);
            }
          gsl_vector_set(w->work, iu, 0.0);
          gsl_vector_set(w->work2, iu - 1, 0.0);

          /* construct right hand side */
          for (k = 0; k < iu - 1; ++k)
            {
              gsl_vector_set(w->work, k,
                             -gsl_vector_get(w->work, iu - 1) *
                             gsl_matrix_get(T, k, iu - 1));
              gsl_vector_set(w->work2, k,
                             -gsl_vector_get(w->work2, iu) *
                             gsl_matrix_get(T, k, iu));
            }

          /*
           * We must solve the upper quasi-triangular system:
           *
           * [ T(1:i-2,1:i-2) - lambda*I ] * X = s*(work + i*work2)
           */

          for (l = i - 2; l >= 0; --l)
            {
              size_t lu = (size_t) l;

              if (lu == 0)
                complex_pair = 0;
              else
                complex_pair = gsl_matrix_get(T, lu, lu - 1) != 0.0;

              if (!complex_pair)
                {
                  gsl_complex bval;
                  gsl_complex x;

                  /*
                   * 1-by-1 diagonal block - solve the system:
                   *
                   * (T_{ll} - lambda)*x = work + i*work2
                   */

                  Tv = gsl_matrix_submatrix(T, lu, lu, 1, 1);
                  bv = gsl_vector_complex_view_array(dat, 1);
                  xv = gsl_vector_complex_view_array(dat_X, 1);

                  GSL_SET_COMPLEX(&bval,
                                  gsl_vector_get(w->work, lu),
                                  gsl_vector_get(w->work2, lu));
                  gsl_vector_complex_set(&bv.vector, 0, bval);

                  nonsymmv_solve_equation_z(&Tv.matrix,
                                            &lambda,
                                            &bv.vector,
                                            &xv.vector,
                                            &scale,
                                            &xnorm,
                                            smin);

                  if (xnorm > 1.0)
                    {
                      if (gsl_vector_get(w->work3, lu) > bignum / xnorm)
                        {
                          gsl_blas_zdscal(1.0/xnorm, &xv.vector);
                          scale /= xnorm;
                        }
                    }

                  /* scale if necessary */
                  if (scale != 1.0)
                    {
                      gsl_vector_view wv;

                      wv = gsl_vector_subvector(w->work, 0, iu + 1);
                      gsl_blas_dscal(scale, &wv.vector);
                      wv = gsl_vector_subvector(w->work2, 0, iu + 1);
                      gsl_blas_dscal(scale, &wv.vector);
                    }

                  x = gsl_vector_complex_get(&xv.vector, 0);
                  gsl_vector_set(w->work, lu, GSL_REAL(x));
                  gsl_vector_set(w->work2, lu, GSL_IMAG(x));

                  /* update the right hand side */
                  if (lu > 0)
                    {
                      gsl_vector_view v1, v2;

                      v1 = gsl_matrix_column(T, lu);
                      v1 = gsl_vector_subvector(&v1.vector, 0, lu);
                      v2 = gsl_vector_subvector(w->work, 0, lu);
                      gsl_blas_daxpy(-GSL_REAL(x), &v1.vector, &v2.vector);

                      v2 = gsl_vector_subvector(w->work2, 0, lu);
                      gsl_blas_daxpy(-GSL_IMAG(x), &v1.vector, &v2.vector);
                    } /* if (lu > 0) */
                } /* if (!complex_pair) */
              else
                {
                  gsl_complex b1, b2, x1, x2;

                  /*
                   * 2-by-2 diagonal block - solve the system
                   */

                  Tv = gsl_matrix_submatrix(T, lu - 1, lu - 1, 2, 2);
                  bv = gsl_vector_complex_view_array(dat, 2);
                  xv = gsl_vector_complex_view_array(dat_X, 2);

                  GSL_SET_COMPLEX(&b1,
                                  gsl_vector_get(w->work, lu - 1),
                                  gsl_vector_get(w->work2, lu - 1));
                  GSL_SET_COMPLEX(&b2,
                                  gsl_vector_get(w->work, lu),
                                  gsl_vector_get(w->work2, lu));
                  gsl_vector_complex_set(&bv.vector, 0, b1);
                  gsl_vector_complex_set(&bv.vector, 1, b2);

                  nonsymmv_solve_equation_z(&Tv.matrix,
                                            &lambda,
                                            &bv.vector,
                                            &xv.vector,
                                            &scale,
                                            &xnorm,
                                            smin);

                  x1 = gsl_vector_complex_get(&xv.vector, 0);
                  x2 = gsl_vector_complex_get(&xv.vector, 1);

                  if (xnorm > 1.0)
                    {
                      double beta;

                      beta = GSL_MAX(gsl_vector_get(w->work3, lu - 1),
                                     gsl_vector_get(w->work3, lu));
                      if (beta > bignum / xnorm)
                        {
                          gsl_blas_zdscal(1.0/xnorm, &xv.vector);
                          scale /= xnorm;
                        }
                    }

                  /* scale if necessary */
                  if (scale != 1.0)
                    {
                      gsl_vector_view wv;

                      wv = gsl_vector_subvector(w->work, 0, iu + 1);
                      gsl_blas_dscal(scale, &wv.vector);
                      wv = gsl_vector_subvector(w->work2, 0, iu + 1);
                      gsl_blas_dscal(scale, &wv.vector);
                    }
                  gsl_vector_set(w->work, lu - 1, GSL_REAL(x1));
                  gsl_vector_set(w->work, lu, GSL_REAL(x2));
                  gsl_vector_set(w->work2, lu - 1, GSL_IMAG(x1));
                  gsl_vector_set(w->work2, lu, GSL_IMAG(x2));

                  /* update right hand side */
                  if (lu > 1)
                    {
                      gsl_vector_view v1, v2, v3, v4;

                      v1 = gsl_matrix_column(T, lu - 1);
                      v1 = gsl_vector_subvector(&v1.vector, 0, lu - 1);
                      v4 = gsl_matrix_column(T, lu);
                      v4 = gsl_vector_subvector(&v4.vector, 0, lu - 1);
                      v2 = gsl_vector_subvector(w->work, 0, lu - 1);
                      v3 = gsl_vector_subvector(w->work2, 0, lu - 1);

                      gsl_blas_daxpy(-GSL_REAL(x1), &v1.vector, &v2.vector);
                      gsl_blas_daxpy(-GSL_REAL(x2), &v4.vector, &v2.vector);
                      gsl_blas_daxpy(-GSL_IMAG(x1), &v1.vector, &v3.vector);
                      gsl_blas_daxpy(-GSL_IMAG(x2), &v4.vector, &v3.vector);
                    } /* if (lu > 1) */

                  --l;
                } /* if (complex_pair) */
            } /* for (l = i - 2; l >= 0; --l) */

          /*
           * At this point, work + i*work2 is an eigenvector
           * of T - backtransform to get an eigenvector of the
           * original matrix
           */

          y = gsl_matrix_column(Z, iu - 1);
          y2 = gsl_matrix_column(Z, iu);

          if (iu > 1)
            {
              gsl_vector_view x;

              /* compute real part of eigenvectors */

              Zv = gsl_matrix_submatrix(Z, 0, 0, N, iu - 1);
              x = gsl_vector_subvector(w->work, 0, iu - 1);

              gsl_blas_dgemv(CblasNoTrans,
                             1.0,
                             &Zv.matrix,
                             &x.vector,
                             gsl_vector_get(w->work, iu - 1),
                             &y.vector);

              /* now compute the imaginary part */
              x = gsl_vector_subvector(w->work2, 0, iu - 1);

              gsl_blas_dgemv(CblasNoTrans,
                             1.0,
                             &Zv.matrix,
                             &x.vector,
                             gsl_vector_get(w->work2, iu),
                             &y2.vector);
            }
          else
            {
              gsl_blas_dscal(gsl_vector_get(w->work, iu - 1), &y.vector);
              gsl_blas_dscal(gsl_vector_get(w->work2, iu), &y2.vector);
            }

          /*
           * Now store the eigenvectors into evec - the real parts
           * are Z(:,iu - 1) and the imaginary parts are
           * +/- Z(:,iu)
           */

          /* get views of the two eigenvector slots */
          ecol = gsl_matrix_complex_column(evec, iu - 1);
          ecol2 = gsl_matrix_complex_column(evec, iu);

          /*
           * save imaginary part first as it may get overwritten
           * when copying the real part due to our storage scheme
           * in Z/evec
           */
          ev = gsl_vector_complex_imag(&ecol.vector);
          ev2 = gsl_vector_complex_imag(&ecol2.vector);
          scale = 0.0;
          for (ii = 0; ii < N; ++ii)
            {
              double a = gsl_vector_get(&y2.vector, ii);

              scale = GSL_MAX(scale,
                              fabs(a) + fabs(gsl_vector_get(&y.vector, ii)));

              gsl_vector_set(&ev.vector, ii, a);
              gsl_vector_set(&ev2.vector, ii, -a);
            }

          /* now save the real part */
          ev = gsl_vector_complex_real(&ecol.vector);
          ev2 = gsl_vector_complex_real(&ecol2.vector);
          for (ii = 0; ii < N; ++ii)
            {
              double a = gsl_vector_get(&y.vector, ii);

              gsl_vector_set(&ev.vector, ii, a);
              gsl_vector_set(&ev2.vector, ii, a);
            }

          if (scale != 0.0)
            scale = 1.0 / scale;

          /* scale by largest element magnitude */

          gsl_blas_zdscal(scale, &ecol.vector);
          gsl_blas_zdscal(scale, &ecol2.vector);

          /*
           * decrement i since we took care of two eigenvalues at
           * the same time
           */
          --i;
        } /* if (GSL_IMAG(lambda) != 0.0) */
    } /* for (i = (int) N - 1; i >= 0; --i) */
} /* nonsymmv_get_right_eigenvectors() */

/*
nonsymmv_solve_equation()

  Solve the equation which comes up in the back substitution
when computing eigenvectors corresponding to real eigenvalues.
The equation that is solved is:

(A - z*I)*x = s*b

where

A is n-by-n with n = 1 or 2
b and x are n-by-1 real vectors
s is a scaling factor set by this function to prevent overflow in x

Inputs: A     - square matrix (n-by-n)
        z     - real scalar (eigenvalue)
        b     - right hand side vector
        x     - (output) where to store solution
        s     - (output) scale factor
        xnorm - (output) infinity norm of X
        smin  - lower bound on singular values of A - if A - z*I
                is less than this value, we'll use smin*I instead.
                This value should be a safe distance above underflow.

Notes: 1) A and b are not changed on output
       2) Based on lapack routine DLALN2
*/

static inline void
nonsymmv_solve_equation(gsl_matrix *A, double z, gsl_vector *b,
                        gsl_vector *x, double *s, double *xnorm,
                        double smin)
{
  size_t N = A->size1;
  double bnorm;
  double scale = 1.0;

  if (N == 1)
    {
      double c,     /* denominator */
             cnorm; /* |c| */

      /*
       * we have a 1-by-1 (real) scalar system to solve:
       *
       * (a - z)*x = b
       * with z real
       */

      /* c = a - z */
      c = gsl_matrix_get(A, 0, 0) - z;
      cnorm = fabs(c);

      if (cnorm < smin)
        {
          /* set c = smin*I */
          c = smin;
          cnorm = smin;
        }

      /* check scaling for x = b / c */
      bnorm = fabs(gsl_vector_get(b, 0));
      if (cnorm < 1.0 && bnorm > 1.0)
        {
          if (bnorm > GSL_NONSYMMV_BIGNUM*cnorm)
            scale = 1.0 / bnorm;
        }

      /* compute x */
      gsl_vector_set(x, 0, gsl_vector_get(b, 0) * scale / c);
      *xnorm = fabs(gsl_vector_get(x, 0));
    } /* if (N == 1) */
  else
    {
      double cr[2][2];
      double *crv;
      double cmax;
      size_t icmax, j;
      double bval1, bval2;
      double ur11, ur12, ur22, ur11r;
      double cr21, cr22;
      double lr21;
      double b1, b2, bbnd;
      double x1, x2;
      double temp;
      size_t ipivot[4][4] = { { 0, 1, 2, 3 },
                              { 1, 0, 3, 2 },
                              { 2, 3, 0, 1 },
                              { 3, 2, 1, 0 } };
      int rswap[4] = { 0, 1, 0, 1 };
      int zswap[4] = { 0, 0, 1, 1 };

      /*
       * we have a 2-by-2 real system to solve:
       *
       * [ A11 - z   A12   ] [ x1 ] = [ b1 ]
       * [   A21   A22 - z ] [ x2 ]   [ b2 ]
       *
       * (z real)
       */

      crv = (double *) cr;

      /*
       * compute the real part of C = A - z*I - use column ordering
       * here since porting from lapack
       */
      cr[0][0] = gsl_matrix_get(A, 0, 0) - z;
      cr[1][1] = gsl_matrix_get(A, 1, 1) - z;
      cr[0][1] = gsl_matrix_get(A, 1, 0);
      cr[1][0] = gsl_matrix_get(A, 0, 1);

      /* find the largest element in C */
      cmax = 0.0;
      icmax = 0;
      for (j = 0; j < 4; ++j)
        {
          if (fabs(crv[j]) > cmax)
            {
              cmax = fabs(crv[j]);
              icmax = j;
            }
        }

      bval1 = gsl_vector_get(b, 0);
      bval2 = gsl_vector_get(b, 1);

      /* if norm(C) < smin, use smin*I */

      if (cmax < smin)
        {
          bnorm = GSL_MAX(fabs(bval1), fabs(bval2));
          if (smin < 1.0 && bnorm > 1.0)
            {
              if (bnorm > GSL_NONSYMMV_BIGNUM*smin)
                scale = 1.0 / bnorm;
            }
          temp = scale / smin;
          gsl_vector_set(x, 0, temp * bval1);
          gsl_vector_set(x, 1, temp * bval2);
          *xnorm = temp * bnorm;
          *s = scale;
          return;
        }

      /* gaussian elimination with complete pivoting */
      ur11 = crv[icmax];
      cr21 = crv[ipivot[1][icmax]];
      ur12 = crv[ipivot[2][icmax]];
      cr22 = crv[ipivot[3][icmax]];
      ur11r = 1.0 / ur11;
      lr21 = ur11r * cr21;
      ur22 = cr22 - ur12 * lr21;

      /* if smaller pivot < smin, use smin */
      if (fabs(ur22) < smin)
        ur22 = smin;

      if (rswap[icmax])
        {
          b1 = bval2;
          b2 = bval1;
        }
      else
        {
          b1 = bval1;
          b2 = bval2;
        }

      b2 -= lr21 * b1;
      bbnd = GSL_MAX(fabs(b1 * (ur22 * ur11r)), fabs(b2));
      if (bbnd > 1.0 && fabs(ur22) < 1.0)
        {
          if (bbnd >= GSL_NONSYMMV_BIGNUM * fabs(ur22))
            scale = 1.0 / bbnd;
        }

      x2 = (b2 * scale) / ur22;
      x1 = (scale * b1) * ur11r - x2 * (ur11r * ur12);
      if (zswap[icmax])
        {
          gsl_vector_set(x, 0, x2);
          gsl_vector_set(x, 1, x1);
        }
      else
        {
          gsl_vector_set(x, 0, x1);
          gsl_vector_set(x, 1, x2);
        }

      *xnorm = GSL_MAX(fabs(x1), fabs(x2));

      /* further scaling if norm(A) norm(X) > overflow */
      if (*xnorm > 1.0 && cmax > 1.0)
        {
          if (*xnorm > GSL_NONSYMMV_BIGNUM / cmax)
            {
              temp = cmax / GSL_NONSYMMV_BIGNUM;
              gsl_blas_dscal(temp, x);
              *xnorm *= temp;
              scale *= temp;
            }
        }
    } /* if (N == 2) */

  *s = scale;
} /* nonsymmv_solve_equation() */

/*
nonsymmv_solve_equation_z()

  Solve the equation which comes up in the back substitution
when computing eigenvectors corresponding to complex eigenvalues.
The equation that is solved is:

(A - z*I)*x = s*b

where

A is n-by-n with n = 1 or 2
b and x are n-by-1 complex vectors
s is a scaling factor set by this function to prevent overflow in x

Inputs: A     - square matrix (n-by-n)
        z     - complex scalar (eigenvalue)
        b     - right hand side vector
        x     - (output) where to store solution
        s     - (output) scale factor
        xnorm - (output) infinity norm of X
        smin  - lower bound on singular values of A - if A - z*I
                is less than this value, we'll use smin*I instead.
                This value should be a safe distance above underflow.

Notes: 1) A and b are not changed on output
       2) Based on lapack routine DLALN2
*/

static inline void
nonsymmv_solve_equation_z(gsl_matrix *A, gsl_complex *z,
                          gsl_vector_complex *b, gsl_vector_complex *x,
                          double *s, double *xnorm, double smin)
{
  size_t N = A->size1;
  double scale = 1.0;
  double bnorm;

  if (N == 1)
    {
      double cr,    /* denominator */
             ci,
             cnorm; /* |c| */
      gsl_complex bval, c, xval, tmp;

      /*
       * we have a 1-by-1 (complex) scalar system to solve:
       *
       * (a - z)*x = b
       * (z is complex, a is real)
       */

      /* c = a - z */
      cr = gsl_matrix_get(A, 0, 0) - GSL_REAL(*z);
      ci = -GSL_IMAG(*z);
      cnorm = fabs(cr) + fabs(ci);

      if (cnorm < smin)
        {
          /* set c = smin*I */
          cr = smin;
          ci = 0.0;
          cnorm = smin;
        }

      /* check scaling for x = b / c */
      bval = gsl_vector_complex_get(b, 0);
      bnorm = fabs(GSL_REAL(bval)) + fabs(GSL_IMAG(bval));
      if (cnorm < 1.0 && bnorm > 1.0)
        {
          if (bnorm > GSL_NONSYMMV_BIGNUM*cnorm)
            scale = 1.0 / bnorm;
        }

      /* compute x */
      GSL_SET_COMPLEX(&tmp, scale*GSL_REAL(bval), scale*GSL_IMAG(bval));
      GSL_SET_COMPLEX(&c, cr, ci);
      xval = gsl_complex_div(tmp, c);

      gsl_vector_complex_set(x, 0, xval);

      *xnorm = fabs(GSL_REAL(xval)) + fabs(GSL_IMAG(xval));
    } /* if (N == 1) */
  else
    {
      double cr[2][2], ci[2][2];
      double *civ, *crv;
      double cmax;
      gsl_complex bval1, bval2;
      gsl_complex xval1, xval2;
      double xr1, xi1;
      size_t icmax;
      size_t j;
      double temp;
      double ur11, ur12, ur22, ui11, ui12, ui22, ur11r, ui11r;
      double ur12s, ui12s;
      double u22abs;
      double lr21, li21;
      double cr21, cr22, ci21, ci22;
      double br1, bi1, br2, bi2, bbnd;
      gsl_complex b1, b2;
      size_t ipivot[4][4] = { { 0, 1, 2, 3 },
                              { 1, 0, 3, 2 },
                              { 2, 3, 0, 1 },
                              { 3, 2, 1, 0 } };
      int rswap[4] = { 0, 1, 0, 1 };
      int zswap[4] = { 0, 0, 1, 1 };

      /*
       * complex 2-by-2 system:
       *
       * [ A11 - z   A12   ] [ X1 ] = [ B1 ]
       * [   A21   A22 - z ] [ X2 ]   [ B2 ]
       *
       * (z complex)
       *
       * where the X and B values are complex.
       */

      civ = (double *) ci;
      crv = (double *) cr;

      /*
       * compute the real part of C = A - z*I - use column ordering
       * here since porting from lapack
       */
      cr[0][0] = gsl_matrix_get(A, 0, 0) - GSL_REAL(*z);
      cr[1][1] = gsl_matrix_get(A, 1, 1) - GSL_REAL(*z);
      cr[0][1] = gsl_matrix_get(A, 1, 0);
      cr[1][0] = gsl_matrix_get(A, 0, 1);

      /* compute the imaginary part */
      ci[0][0] = -GSL_IMAG(*z);
      ci[0][1] = 0.0;
      ci[1][0] = 0.0;
      ci[1][1] = -GSL_IMAG(*z);

      cmax = 0.0;
      icmax = 0;

      for (j = 0; j < 4; ++j)
        {
          if (fabs(crv[j]) + fabs(civ[j]) > cmax)
            {
              cmax = fabs(crv[j]) + fabs(civ[j]);
              icmax = j;
            }
        }

      bval1 = gsl_vector_complex_get(b, 0);
      bval2 = gsl_vector_complex_get(b, 1);

      /* if norm(C) < smin, use smin*I */
      if (cmax < smin)
        {
          bnorm = GSL_MAX(fabs(GSL_REAL(bval1)) + fabs(GSL_IMAG(bval1)),
                          fabs(GSL_REAL(bval2)) + fabs(GSL_IMAG(bval2)));
          if (smin < 1.0 && bnorm > 1.0)
            {
              if (bnorm > GSL_NONSYMMV_BIGNUM*smin)
                scale = 1.0 / bnorm;
            }

          temp = scale / smin;
          xval1 = gsl_complex_mul_real(bval1, temp);
          xval2 = gsl_complex_mul_real(bval2, temp);
          gsl_vector_complex_set(x, 0, xval1);
          gsl_vector_complex_set(x, 1, xval2);
          *xnorm = temp * bnorm;
          *s = scale;
          return;
        }

      /* gaussian elimination with complete pivoting */
      ur11 = crv[icmax];
      ui11 = civ[icmax];
      cr21 = crv[ipivot[1][icmax]];
      ci21 = civ[ipivot[1][icmax]];
      ur12 = crv[ipivot[2][icmax]];
      ui12 = civ[ipivot[2][icmax]];
      cr22 = crv[ipivot[3][icmax]];
      ci22 = civ[ipivot[3][icmax]];

      if (icmax == 0 || icmax == 3)
        {
          /* off diagonals of pivoted C are real */
          if (fabs(ur11) > fabs(ui11))
            {
              temp = ui11 / ur11;
              ur11r = 1.0 / (ur11 * (1.0 + temp*temp));
              ui11r = -temp * ur11r;
            }
          else
            {
              temp = ur11 / ui11;
              ui11r = -1.0 / (ui11 * (1.0 + temp*temp));
              ur11r = -temp*ui11r;
            }
          lr21 = cr21 * ur11r;
          li21 = cr21 * ui11r;
          ur12s = ur12 * ur11r;
          ui12s = ur12 * ui11r;
          ur22 = cr22 - ur12 * lr21;
          ui22 = ci22 - ur12 * li21;
        }
      else
        {
          /* diagonals of pivoted C are real */
          ur11r = 1.0 / ur11;
          ui11r = 0.0;
          lr21 = cr21 * ur11r;
          li21 = ci21 * ur11r;
          ur12s = ur12 * ur11r;
          ui12s = ui12 * ur11r;
          ur22 = cr22 - ur12 * lr21 + ui12 * li21;
          ui22 = -ur12 * li21 - ui12 * lr21;
        }

      u22abs = fabs(ur22) + fabs(ui22);

      /* if smaller pivot < smin, use smin */
      if (u22abs < smin)
        {
          ur22 = smin;
          ui22 = 0.0;
        }

      if (rswap[icmax])
        {
          br2 = GSL_REAL(bval1);
          bi2 = GSL_IMAG(bval1);
          br1 = GSL_REAL(bval2);
          bi1 = GSL_IMAG(bval2);
        }
      else
        {
          br1 = GSL_REAL(bval1);
          bi1 = GSL_IMAG(bval1);
          br2 = GSL_REAL(bval2);
          bi2 = GSL_IMAG(bval2);
        }

      br2 += li21*bi1 - lr21*br1;
      bi2 -= li21*br1 + lr21*bi1;
      bbnd = GSL_MAX((fabs(br1) + fabs(bi1)) *
                     (u22abs * (fabs(ur11r) + fabs(ui11r))),
                     fabs(br2) + fabs(bi2));
      if (bbnd > 1.0 && u22abs < 1.0)
        {
          if (bbnd >= GSL_NONSYMMV_BIGNUM*u22abs)
            {
              scale = 1.0 / bbnd;
              br1 *= scale;
              bi1 *= scale;
              br2 *= scale;
              bi2 *= scale;
            }
        }

      GSL_SET_COMPLEX(&b1, br2, bi2);
      GSL_SET_COMPLEX(&b2, ur22, ui22);
      xval2 = gsl_complex_div(b1, b2);

      xr1 = ur11r*br1 - ui11r*bi1 - ur12s*GSL_REAL(xval2) + ui12s*GSL_IMAG(xval2);
      xi1 = ui11r*br1 + ur11r*bi1 - ui12s*GSL_REAL(xval2) - ur12s*GSL_IMAG(xval2);
      GSL_SET_COMPLEX(&xval1, xr1, xi1);

      if (zswap[icmax])
        {
          gsl_vector_complex_set(x, 0, xval2);
          gsl_vector_complex_set(x, 1, xval1);
        }
      else
        {
          gsl_vector_complex_set(x, 0, xval1);
          gsl_vector_complex_set(x, 1, xval2);
        }

      *xnorm = GSL_MAX(fabs(GSL_REAL(xval1)) + fabs(GSL_IMAG(xval1)),
                       fabs(GSL_REAL(xval2)) + fabs(GSL_IMAG(xval2)));

      /* further scaling if norm(A) norm(X) > overflow */
      if (*xnorm > 1.0 && cmax > 1.0)
        {
          if (*xnorm > GSL_NONSYMMV_BIGNUM / cmax)
            {
              temp = cmax / GSL_NONSYMMV_BIGNUM;
              gsl_blas_zdscal(temp, x);
              *xnorm *= temp;
              scale *= temp;
            }
        }
    } /* if (N == 2) */

  *s = scale;
} /* nonsymmv_solve_equation_z() */

/*
nonsymmv_normalize_eigenvectors()
  Normalize eigenvectors so that their Euclidean norm is 1

Inputs: eval - eigenvalues
        evec - eigenvectors
*/

static void
nonsymmv_normalize_eigenvectors(gsl_vector_complex *eval,
                                gsl_matrix_complex *evec)
{
  const size_t N = evec->size1;
  size_t i;     /* looping */
  gsl_complex ei;
  gsl_vector_complex_view vi;
  gsl_vector_view re, im;
  double scale; /* scaling factor */

  for (i = 0; i < N; ++i)
    {
      ei = gsl_vector_complex_get(eval, i);
      vi = gsl_matrix_complex_column(evec, i);

      re = gsl_vector_complex_real(&vi.vector);

      if (GSL_IMAG(ei) == 0.0)
        {
          scale = 1.0 / gsl_blas_dnrm2(&re.vector);
          gsl_blas_dscal(scale, &re.vector);
        }
      else if (GSL_IMAG(ei) > 0.0)
        {
          im = gsl_vector_complex_imag(&vi.vector);

          scale = 1.0 / gsl_hypot(gsl_blas_dnrm2(&re.vector),
                                  gsl_blas_dnrm2(&im.vector));
          gsl_blas_zdscal(scale, &vi.vector);

          vi = gsl_matrix_complex_column(evec, i + 1);
          gsl_blas_zdscal(scale, &vi.vector);
        }
    }
} /* nonsymmv_normalize_eigenvectors() */

max_stars_repo_name: manggoguy/parsec-modified

max_issues_repo_name: manggoguy/parsec-modified

max_forks_repo_name: manggoguy/parsec-modified


text:
/////////////////////////////////////////////////////////////////////
// = NMatrix
//
// A linear algebra library for scientific computation in Ruby.
// NMatrix is part of SciRuby.
//
// NMatrix was originally inspired by and derived from NArray, by
// Masahiro Tanaka: http://narray.rubyforge.org
//
// == Copyright Information
//
// SciRuby is Copyright (c) 2010 - 2014, Ruby Science Foundation
// NMatrix is Copyright (c) 2012 - 2014, John Woods and the Ruby Science Foundation
//
// Please see LICENSE.txt for additional copyright notices.
//
// == Contributing
//
// By contributing source code to SciRuby, you agree to be bound by
// our Contributor Agreement:
//
// * https://github.com/SciRuby/sciruby/wiki/Contributor-Agreement
//
// == gemm.h
//
// Header file for interface with ATLAS's CBLAS gemm functions and
// native templated version of LAPACK's gemm function.
//

#ifndef GEMM_H
# define GEMM_H

#if !(defined(HAVE_CBLAS_H) || defined(HAVE_ATLAS_CBLAS_H))
#include "cblas_enums.h"
#endif

extern "C" { // These need to be in an extern "C" block or you'll get all kinds of undefined symbol errors.
#if defined HAVE_CBLAS_H
  #include <cblas.h>
#elif defined HAVE_ATLAS_CBLAS_H
  #include <atlas/cblas.h>
#endif
}

namespace nm { namespace math {
/*
 * GEneral Matrix Multiplication: based on dgemm.f from Netlib.
 *
 * This is an extremely inefficient algorithm. Recommend using ATLAS' version instead.
 *
 * Template parameters: LT -- long version of type T. Type T is the matrix dtype.
 *
 * This version throws no errors. Use gemm<DType> instead for error checking.
 */
template <typename DType>
inline void gemm_nothrow(const enum CBLAS_TRANSPOSE TransA, const enum CBLAS_TRANSPOSE TransB, const int M, const int N, const int K,
                 const DType* alpha, const DType* A, const int lda, const DType* B, const int ldb, const DType* beta, DType* C, const int ldc)
{

  typename LongDType<DType>::type temp;

  // Quick return if possible
  if (!M or !N or ((*alpha == 0 or !K) and *beta == 1)) return;

  // For alpha = 0
  if (*alpha == 0) {
    if (*beta == 0) {
      for (int j = 0; j < N; ++j)
        for (int i = 0; i < M; ++i) {
          C[i+j*ldc] = 0;
        }
    } else {
      for (int j = 0; j < N; ++j)
        for (int i = 0; i < M; ++i) {
          C[i+j*ldc] *= *beta;
        }
    }
    return;
  }

  // Start the operations
  if (TransB == CblasNoTrans) {
    if (TransA == CblasNoTrans) {
      // C = alpha*A*B+beta*C
      for (int j = 0; j < N; ++j) {
        if (*beta == 0) {
          for (int i = 0; i < M; ++i) {
            C[i+j*ldc] = 0;
          }
        } else if (*beta != 1) {
          for (int i = 0; i < M; ++i) {
            C[i+j*ldc] *= *beta;
          }
        }

        for (int l = 0; l < K; ++l) {
          if (B[l+j*ldb] != 0) {
            temp = *alpha * B[l+j*ldb];
            for (int i = 0; i < M; ++i) {
              C[i+j*ldc] += A[i+l*lda] * temp;
            }
          }
        }
      }

    } else {

      // C = alpha*A**DType*B + beta*C
      for (int j = 0; j < N; ++j) {
        for (int i = 0; i < M; ++i) {
          temp = 0;
          for (int l = 0; l < K; ++l) {
            temp += A[l+i*lda] * B[l+j*ldb];
          }

          if (*beta == 0) {
            C[i+j*ldc] = *alpha*temp;
          } else {
            C[i+j*ldc] = *alpha*temp + *beta*C[i+j*ldc];
          }
        }
      }

    }

  } else if (TransA == CblasNoTrans) {

    // C = alpha*A*B**T + beta*C
    for (int j = 0; j < N; ++j) {
      if (*beta == 0) {
        for (int i = 0; i < M; ++i) {
          C[i+j*ldc] = 0;
        }
      } else if (*beta != 1) {
        for (int i = 0; i < M; ++i) {
          C[i+j*ldc] *= *beta;
        }
      }

      for (int l = 0; l < K; ++l) {
        if (B[j+l*ldb] != 0) {
          temp = *alpha * B[j+l*ldb];
          for (int i = 0; i < M; ++i) {
            C[i+j*ldc] += A[i+l*lda] * temp;
          }
        }
      }

    }

  } else {

    // C = alpha*A**DType*B**T + beta*C
    for (int j = 0; j < N; ++j) {
      for (int i = 0; i < M; ++i) {
        temp = 0;
        for (int l = 0; l < K; ++l) {
          temp += A[l+i*lda] * B[j+l*ldb];
        }

        if (*beta == 0) {
          C[i+j*ldc] = *alpha*temp;
        } else {
          C[i+j*ldc] = *alpha*temp + *beta*C[i+j*ldc];
        }
      }
    }

  }

  return;
}

template <typename DType>
inline void gemm(const enum CBLAS_ORDER Order, const enum CBLAS_TRANSPOSE TransA, const enum CBLAS_TRANSPOSE TransB, const int M, const int N, const int K,
                 const DType* alpha, const DType* A, const int lda, const DType* B, const int ldb, const DType* beta, DType* C, const int ldc)
{
  if (Order == CblasRowMajor) {
    if (TransA == CblasNoTrans) {
      if (lda < std::max(K,1)) {
        rb_raise(rb_eArgError, "lda must be >= MAX(K,1): lda=%d K=%d", lda, K);
      }
    } else {
      if (lda < std::max(M,1)) { // && TransA == CblasTrans
        rb_raise(rb_eArgError, "lda must be >= MAX(M,1): lda=%d M=%d", lda, M);
      }
    }

    if (TransB == CblasNoTrans) {
      if (ldb < std::max(N,1)) {
        rb_raise(rb_eArgError, "ldb must be >= MAX(N,1): ldb=%d N=%d", ldb, N);
      }
    } else {
      if (ldb < std::max(K,1)) {
        rb_raise(rb_eArgError, "ldb must be >= MAX(K,1): ldb=%d K=%d", ldb, K);
      }
    }

    if (ldc < std::max(N,1)) {
      rb_raise(rb_eArgError, "ldc must be >= MAX(N,1): ldc=%d N=%d", ldc, N);
    }
  } else { // CblasColMajor
    if (TransA == CblasNoTrans) {
      if (lda < std::max(M,1)) {
        rb_raise(rb_eArgError, "lda must be >= MAX(M,1): lda=%d M=%d", lda, M);
      }
    } else {
      if (lda < std::max(K,1)) { // && TransA == CblasTrans
        rb_raise(rb_eArgError, "lda must be >= MAX(K,1): lda=%d K=%d", lda, K);
      }
    }

    if (TransB == CblasNoTrans) {
      if (ldb < std::max(K,1)) {
        rb_raise(rb_eArgError, "ldb must be >= MAX(K,1): ldb=%d N=%d", ldb, K);
      }
    } else {
      if (ldb < std::max(N,1)) { // NOTE: This error message is actually wrong in the ATLAS source currently. Or are we wrong?
        rb_raise(rb_eArgError, "ldb must be >= MAX(N,1): ldb=%d N=%d", ldb, N);
      }
    }

    if (ldc < std::max(M,1)) {
      rb_raise(rb_eArgError, "ldc must be >= MAX(M,1): ldc=%d N=%d", ldc, M);
    }
  }

  /*
   * Call SYRK when that's what the user is actually asking for; just handle beta=0, because beta=X requires
   * we copy C and then subtract to preserve asymmetry.
   */

  if (A == B && M == N && TransA != TransB && lda == ldb && beta == 0) {
    rb_raise(rb_eNotImpError, "syrk and syreflect not implemented");
    /*syrk<DType>(CblasUpper, (Order == CblasColMajor) ? TransA : TransB, N, K, alpha, A, lda, beta, C, ldc);
    syreflect(CblasUpper, N, C, ldc);
    */
  }

  if (Order == CblasRowMajor)    gemm_nothrow<DType>(TransB, TransA, N, M, K, alpha, B, ldb, A, lda, beta, C, ldc);
  else                           gemm_nothrow<DType>(TransA, TransB, M, N, K, alpha, A, lda, B, ldb, beta, C, ldc);

}

#if (defined(HAVE_CBLAS_H) || defined(HAVE_ATLAS_CBLAS_H))
template <>
inline void gemm(const enum CBLAS_ORDER Order, const enum CBLAS_TRANSPOSE TransA, const enum CBLAS_TRANSPOSE TransB, const int M, const int N, const int K,
          const float* alpha, const float* A, const int lda, const float* B, const int ldb, const float* beta, float* C, const int ldc) {
  cblas_sgemm(Order, TransA, TransB, M, N, K, *alpha, A, lda, B, ldb, *beta, C, ldc);
}

template <>
inline void gemm(const enum CBLAS_ORDER Order, const enum CBLAS_TRANSPOSE TransA, const enum CBLAS_TRANSPOSE TransB, const int M, const int N, const int K,
          const double* alpha, const double* A, const int lda, const double* B, const int ldb, const double* beta, double* C, const int ldc) {
  cblas_dgemm(Order, TransA, TransB, M, N, K, *alpha, A, lda, B, ldb, *beta, C, ldc);
}

template <>
inline void gemm(const enum CBLAS_ORDER Order, const enum CBLAS_TRANSPOSE TransA, const enum CBLAS_TRANSPOSE TransB, const int M, const int N, const int K,
          const Complex64* alpha, const Complex64* A, const int lda, const Complex64* B, const int ldb, const Complex64* beta, Complex64* C, const int ldc) {
  cblas_cgemm(Order, TransA, TransB, M, N, K, alpha, A, lda, B, ldb, beta, C, ldc);
}

template <>
inline void gemm(const enum CBLAS_ORDER Order, const enum CBLAS_TRANSPOSE TransA, const enum CBLAS_TRANSPOSE TransB, const int M, const int N, const int K,
          const Complex128* alpha, const Complex128* A, const int lda, const Complex128* B, const int ldb, const Complex128* beta, Complex128* C, const int ldc) {
  cblas_zgemm(Order, TransA, TransB, M, N, K, alpha, A, lda, B, ldb, beta, C, ldc);
}
#endif

}} // end of namespace nm::math

#endif // GEMM_H

max_stars_repo_name: cjfuller/nmatrix

max_issues_repo_name: cjfuller/nmatrix

max_forks_repo_name: cjfuller/nmatrix


text:
/* min/test_funcs.c
 *
 * Copyright (C) 1996, 1997, 1998, 1999, 2000 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 2 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 <gsl/gsl_math.h>
#include <gsl/gsl_min.h>

#include "test.h"

gsl_function create_function (double (*f)(double, void *))
{
  gsl_function F ;
  F.function = f ;
  F.params = 0 ;
  return F ;
}

double
f_cos (double x, void * p)
{
  p = 0;  /* avoid warning about unused parameter */
  return cos(x);
}

/* f(x) = x^4 - 1 */
/* minimum at x = 0 */

double
func1 (double x, void * p)
{
  p = 0;  /* avoid warning about unused parameter */
  return pow (x, 4.0) - 1;
}

/* f(x) = sqrt(|x|) */
/* minimum at x = 0 */

double
func2 (double x, void * p)
{
  p = 0;  /* avoid warning about unused parameter */
  return sqrt(fabs(x));
}

/* f(x) = 1 for x < 1 and -exp(-x) for x >= 1 */
/* minimum at x = 1 */

double
func3 (double x, void * p)
{
  p = 0;  /* avoid warning about unused parameter */

  if (x < 1)
    return 1 ;
  else
    return - exp(-x) ;
}

/* f(x) = x - 30/(1+1e5*(x-0.8)**2) */
/* minimum near x = 0.8 */

double
func4 (double x, void * p)
{
  p = 0;  /* avoid warning about unused parameter */

  return x - 30.0 / (1.0 + 1e5 * pow(x-0.8, 2.0));
}