text:
using Test, ModelConstructors, DSGE, Dates, FileIO, Random

generate_fulldist_forecast_data = false
if VERSION < v"1.5"
    ver = "111"
else
    ver = "150"
end

Random.seed!(1793)

@testset "Test regime switching" begin
    n_reg_temp = 14

    m = Model1002("ss10", custom_settings = Dict{Symbol, Setting}(:add_pgap => Setting(:add_pgap, true)))
    m <= Setting(:date_forecast_start, Date(2020, 6, 30))
    m <= Setting(:date_conditional_end, Date(2020, 6, 30))
    m <= Setting(:regime_switching, true)
    m <= Setting(:regime_dates, Dict{Int, Date}(1 => date_presample_start(m),
                                                  2 => Date(2020, 3, 31),
                                                  3 => Date(2020, 6, 30)))
    m = setup_regime_switching_inds!(m)

    t_s, r_s, c_s = solve(m, regimes = Int[3])

    regime_dates = Dict{Int, Date}()
    regime_dates[1] = date_presample_start(m)
    for (i, date) in zip(2:n_reg_temp, Date(2020,3,31):Dates.Month(3):Date(3000,3,31))
        regime_dates[i] = date
    end
    m <= Setting(:regime_dates, regime_dates)

    m <= Setting(:regime_switching, true)

    m = setup_regime_switching_inds!(m)

    # Use gensy2 for temporary rule
    m <= Setting(:gensys2, true)
    # Replace eqcond with temp rule
    m <= Setting(:replace_eqcond, true)
    # Which rule to replace with inn nwhich periods
    replace_eqcond = Dict{Int, Function}()
    for i in 3:n_reg_temp
       replace_eqcond[i] = ngdp_replace_eq_entries
   end

    m <= Setting(:replace_eqcond_func_dict, replace_eqcond)

    m <= Setting(:pgap_value, 12.0 : 0.0)
    m <= Setting(:pgap_type, :ngdp)
    m = setup_regime_switching_inds!(m)

    sys_temp = solve(m, regime_switching = true, regimes = collect(1:n_reg_temp), hist_regimes = [1, 2],
                     fcast_regimes = collect(3:n_reg_temp))
    Γ0s = Vector{Matrix{Float64}}(undef, n_reg_temp-2)
    Γ1s = Vector{Matrix{Float64}}(undef, n_reg_temp-2)
    Cs = Vector{Vector{Float64}}(undef,  n_reg_temp-2)
    Ψs = Vector{Matrix{Float64}}(undef, n_reg_temp-2)
    Πs = Vector{Matrix{Float64}}(undef, n_reg_temp-2)
    n_endo = length(collect(keys(m.endogenous_states)))
    for i in (3-2):(n_reg_temp-2)
        global Γ0s[i], Γ1s[i], Cs[i], Ψs[i], Πs[i] = eqcond(m, i+2)
    end
    Tcal,Rcal, Ccal = DSGE.gensys_cplus(m, Γ0s, Γ1s, Cs, Ψs, Πs, t_s[1:n_endo, 1:n_endo], r_s[1:n_endo, :], c_s[1:n_endo])
    @test Tcal[1] ≈ sys_temp[1][3][1:n_endo, 1:n_endo]
    # Then when start recrusion with rule, should get same TTTs in every perod
    t, r, c = ngdp_solve(m, regime_switching = false, regimes = Int[3])

    Tcal,Rcal, Ccal = DSGE.gensys_cplus(m, Γ0s, Γ1s, Cs, Ψs, Πs, t[1:n_endo, 1:n_endo], r[1:n_endo, :], c[1:n_endo])
    for i in 1:length(Tcal)
        @test Tcal[i] ≈ t[1:n_endo, 1:n_endo]
    end
end

@testset "Test regime switching IRFs match Forward Guidance method IRFs" begin
    horizon = 60
    m = Model1002()
    m <= Setting(:date_forecast_start, Date(2020, 6, 30))
    m <= Setting(:date_conditional_end, Date(2020, 6, 30))
    m <= Setting(:replace_eqcond, false)
    m <= Setting(:gensys2, false)
    m <= Setting(:regime_switching, false)
    m[:ρ_rm] = 0.0

    system = compute_system(m)
    H = 2
    nstates = size(system[:TTT], 1)
    nshocks = size(system[:RRR], 2)
    nobs         = size(system[:ZZ], 1)
    npseudo      = size(system[:ZZ_pseudo], 1)

    shocks = m.exogenous_shocks

    states = zeros(nstates, horizon) #, nshocks)
    obs    = zeros(nobs,    horizon) #, nshocks)
    pseudo = zeros(npseudo, horizon) #, nshocks)

    s_0 = zeros(nstates)
    ### IRFs to FFR peg for H periods
    # Set anticipated shocks to peg FFR from t+1 to t+H
    FFRpeg = -m[:Rstarn].value #-0.25/4 #/400

    PsiR1 = 0 #constant
    # ZZ matrix
    PsiR2 = zeros(nstates)
    PsiR2[m.endogenous_states[:R_t]] = 1

    Rht = system[:RRR][:,vcat(shocks[:rm_sh], shocks[:rm_shl1]:shocks[Symbol("rm_shl$H")])]

    bb = zeros(H+1,1)
    MH = zeros(H+1,H+1)
    for hh = 1:H+1
        bb[hh,1] = (FFRpeg - PsiR1 - PsiR2'*(system[:TTT])^hh*s_0)
        MH[hh,:] = PsiR2'*(system[:TTT])^(hh-1)*Rht
    end
    monshocks = MH\bb
    #% Simulate the model
    etpeg = zeros(nshocks,horizon)
    etpeg[vcat(shocks[:rm_sh], shocks[:rm_shl1]:shocks[Symbol("rm_shl$H")]),1] = monshocks
    states[:, :], obs[:, :], pseudo[:, :] = forecast(system, s_0, etpeg, enforce_zlb = false)

    m = Model1002()
    m <= Setting(:date_forecast_start, Date(2020, 6, 30))
    m <= Setting(:gensys2, true)
    m <= Setting(:replace_eqcond, true)
    replace_eqcond = Dict{Int, Function}()
    for i in 3:5
        replace_eqcond[i] = zero_rate_replace_eq_entries
    end
    m <= Setting(:replace_eqcond_func_dict, replace_eqcond)

    m <= Setting(:regime_switching, true)
    m <= Setting(:regime_dates, Dict{Int, Date}(1 => date_presample_start(m),
                                                2 => Date(1990, 3, 31),
                                                3 => Date(2020, 6, 30),
                                                4 => Date(2020, 9, 30),
                                                5 => Date(2020, 12, 31),
                                                6 => Date(2021, 3, 31)))

    m = setup_regime_switching_inds!(m)

    sys = compute_system(m)
    @test !isapprox(sys[4, :TTT], system[:TTT], atol = 1e-5) # Matrix in regime 4 should not match the matrix w/out regime-switching

    s, o, p = forecast(m, sys, zeros(84), zeros(24, 60),  enforce_zlb = false)

    @test all(isapprox.(s[1:20, :], states[1:20, :], atol = 1e-3))
    @test all(isapprox.(o[1:9, :], obs[1:9, :], atol = 1e-3))
    @test all(isapprox.(p[1:10, :], pseudo[1:10, :], atol = 1e-3))
end

@testset "Test rule switching forecast" begin
    output_vars = [:forecastobs, :histobs, :histpseudo, :forecastpseud]

    m = Model1002("ss10", custom_settings = Dict{Symbol, Setting}(:add_pgap =>
                                                                  Setting(:add_pgap, true)))
    m <= Setting(:date_forecast_start, Date(2020, 6, 30))

    m <= Setting(:regime_switching, true)
    m <= Setting(:regime_dates, Dict{Int, Date}(1 => date_presample_start(m),
                                                  2 => Date(2020, 3, 31),
                                                  3 => Date(2020, 6, 30)))
    m = setup_regime_switching_inds!(m)

    fp = dirname(@__FILE__)
    df = load(joinpath(fp, "../reference/regime_switch_data.jld2"), "regime_switch_df_none")

    m <= Setting(:replace_eqcond, false)
    m <= Setting(:replace_eqcond_func_dict, Dict{Int, Function}(
        3 => ngdp_replace_eq_entries))
    m <= Setting(:pgap_type, :ngdp)
    m <= Setting(:pgap_value, 12.0)

    m <= Setting(:gensys2, false)

    m <= Setting(:alternative_policy, AltPolicy(:ngdp, ngdp_eqcond, ngdp_solve, forecast_init = DSGE.ngdp_forecast_init))

    fcast_altperm = DSGE.forecast_one_draw(m, :mode, :full, output_vars, map(x -> x.value, m.parameters),
                                           df; regime_switching = true, n_regimes = get_setting(m, :n_regimes))

    # Testing ore basic permanent BGDP
    m = Model1002("ss10"; custom_settings = Dict{Symbol, Setting}(:add_pgap =>
                                                                    Setting(:add_pgap, true)))
    m <= Setting(:date_forecast_start, Date(2020, 6, 30))
    m <= Setting(:regime_switching, true)
    m <= Setting(:regime_dates, Dict{Int, Date}(1 => date_presample_start(m),
                                                2 => Date(2020, 3, 31),
                                                3 => Date(2020, 6, 30)))

    m = setup_regime_switching_inds!(m)

    m <= Setting(:pgap_value, 12.0)
    m <= Setting(:pgap_type, :ngdp)

    m <= Setting(:alternative_policy, AltPolicy(:ngdp, ngdp_eqcond, ngdp_solve, forecast_init = DSGE.ngdp_forecast_init))

    fcast_altperm2 = DSGE.forecast_one_draw(m, :mode, :full, output_vars, map(x -> x.value, m.parameters),
                                            df; regime_switching = true, n_regimes = get_setting(m, :n_regimes))

    @test fcast_altperm[:forecastobs] == fcast_altperm2[:forecastobs]
end

@testset "Temporary alternative policies with non-trivial conditional forecasting" begin
    output_vars = [:forecastobs]

    m = Model1002("ss10", custom_settings = Dict{Symbol, Setting}(:add_pgap =>
                                                                  Setting(:add_pgap, true)))

    m <= Setting(:regime_switching, true)
    m <= Setting(:regime_dates, Dict{Int, Date}(1 => date_presample_start(m),
                                                2 => Date(2020, 3, 31),
                                                3 => Date(2020, 6, 30),
                                                4 => Date(2020, 9, 30),
                                                5 => Date(2020, 12, 31),
                                                6 => Date(2021, 3, 31)))
    m <= Setting(:date_conditional_end, Date(2020, 6, 30))
    m <= Setting(:date_forecast_start, Date(2020, 6, 30))
    m <= Setting(:forecast_horizons, 10)
    setup_regime_switching_inds!(m; cond_type = :full)

    fp = dirname(@__FILE__)
    df = load(joinpath(fp, "../reference/regime_switch_data.jld2"), "regime_switch_df_full")
    date_ind = findfirst(df[!, :date] .== Date(2020, 6, 30))
    df[date_ind, :obs_nominalrate] = .3 / 4.
    m <= Setting(:replace_eqcond, true)
    m <= Setting(:replace_eqcond_func_dict, Dict{Int, Function}(
        3 => zero_rate_replace_eq_entries))
    m <= Setting(:pgap_type, :ngdp)
    m <= Setting(:pgap_value, 12.0)
    m <= Setting(:gensys2, true)

    # Check zero rate rule causes an error if specified in a conditional period
    @test_throws DSGE.GensysError DSGE.forecast_one_draw(m, :mode, :full, output_vars, map(x -> x.value, m.parameters),
                                   df; regime_switching = true, n_regimes = get_setting(m, :n_regimes))

    m <= Setting(:replace_eqcond_func_dict, Dict{Int, Function}(
        4 => zero_rate_replace_eq_entries, 5 => zero_rate_replace_eq_entries))
    fcast = DSGE.forecast_one_draw(m, :mode, :full, output_vars, map(x -> x.value, m.parameters),
                                   df; regime_switching = true, n_regimes = get_setting(m, :n_regimes))
    @test fcast[:forecastobs][m.observables[:obs_nominalrate], 1] > .01
    @test abs(fcast[:forecastobs][m.observables[:obs_nominalrate], 2]) < 1e-14
    @test abs(fcast[:forecastobs][m.observables[:obs_nominalrate], 3]) < 1e-14
    @test all(abs.(fcast[:forecastobs][m.observables[:obs_nominalrate], 4:end]) .> .01)

    m <= Setting(:date_conditional_end, Date(2020, 6, 30))
    m <= Setting(:date_forecast_start, Date(2020, 3, 31))
    fcast = DSGE.forecast_one_draw(m, :mode, :full, output_vars, map(x -> x.value, m.parameters),
                                   df; regime_switching = true, n_regimes = get_setting(m, :n_regimes))
    @test fcast[:forecastobs][m.observables[:obs_nominalrate], 1] > .01
    @test abs(fcast[:forecastobs][m.observables[:obs_nominalrate], 2]) > .01
    @test abs(fcast[:forecastobs][m.observables[:obs_nominalrate], 3]) < 1e-14
    @test abs(fcast[:forecastobs][m.observables[:obs_nominalrate], 4]) < 1e-14
    @test all(abs.(fcast[:forecastobs][m.observables[:obs_nominalrate], 5:end]) .> .01)

    m <= Setting(:date_conditional_end, Date(2020, 3, 31))
    m <= Setting(:date_forecast_start, Date(2020, 3, 31))
    fcast = DSGE.forecast_one_draw(m, :mode, :full, output_vars, map(x -> x.value, m.parameters),
                                   df; regime_switching = true, n_regimes = get_setting(m, :n_regimes))
    @test fcast[:forecastobs][m.observables[:obs_nominalrate], 1] > 0.
    @test abs(fcast[:forecastobs][m.observables[:obs_nominalrate], 2]) < 1e-14
    @test abs(fcast[:forecastobs][m.observables[:obs_nominalrate], 3]) < 1e-14
    @test abs(fcast[:forecastobs][m.observables[:obs_nominalrate], 4]) < 1e-14
    @test all(abs.(fcast[:forecastobs][m.observables[:obs_nominalrate], 5:end]) .> 0.)
end

@testset "Calculation of regime indices" begin
    m = AnSchorfheide("ss0")

    rss = Vector{Dict{Int, Date}}(undef, 0)
    fss = Vector{Date}(undef, 0)
    ces = Vector{Date}(undef, 0)
    answers = Dict()
    answers[:full] = []
    answers[:none] = []

    push!(rss, Dict{Int, Date}(1 => date_presample_start(m), 2 => Date(2020, 3, 31), 3 => Date(2020, 6, 30), 4 => Date(2020, 9, 30)))
    push!(fss, Date(2020, 6, 30))
    push!(ces, Date(2020, 6, 30))
    push!(answers[:none], [1:1, 2:60])
    push!(answers[:full], [1:59])

    push!(rss, Dict{Int, Date}(1 => date_presample_start(m), 2 => Date(2020, 3, 31), 3 => Date(2020, 6, 30),
                               4 => Date(2021, 12, 31)))
    push!(fss, Date(2020, 6, 30))
    push!(ces, Date(2020, 6, 30))
    push!(answers[:none], [1:6, 7:60])
    push!(answers[:full], [1:5, 6:59])

    push!(rss, Dict{Int, Date}(1 => date_presample_start(m), 2 => Date(2020, 3, 31), 3 => Date(2020, 12, 31),
                               4 => Date(2021, 12, 31)))
    push!(fss, Date(2020, 6, 30))
    push!(ces, Date(2020, 12, 31))
    push!(answers[:none], [1:2, 3:6, 7:60])
    push!(answers[:full], [1:3, 4:57])

    push!(rss, Dict{Int, Date}(1 => date_presample_start(m), 2 => Date(2020, 3, 31), 3 => Date(2021, 12, 31),
                           4 => Date(2022, 12, 31)))
    push!(fss, Date(2020, 6, 30))
    push!(ces, Date(2020, 12, 31))
    push!(answers[:none], [1:6, 7:10, 11:60])
    push!(answers[:full], [1:3, 4:7, 8:57])

    push!(rss, Dict{Int, Date}(1 => date_presample_start(m), 2 => Date(2020, 6, 30), 3 => Date(2021, 12, 31),
                               4 => Date(2022, 12, 31)))
    push!(fss, Date(2020, 6, 30))
    push!(ces, Date(2020, 12, 31))
    push!(answers[:none], [1:6, 7:10, 11:60])
    push!(answers[:full], [1:3, 4:7, 8:57])

    for cond_type in [:full, :none]
        for (i, rs, fs, ce) in zip(1:length(rss), rss, fss, ces)
            m <= Setting(:regime_dates, rs)
            m <= Setting(:date_forecast_start, fs)
            m <= Setting(:date_conditional_end, ce)
            setup_regime_switching_inds!(m)
            fcast_reg_inds = DSGE.get_fcast_regime_inds(m, forecast_horizons(m; cond_type = cond_type), cond_type)
            @test fcast_reg_inds == answers[cond_type][i]
        end
    end
end

@testset "Full-distribution forecast with conditioning and regime switching during forecast horizon" begin
    # Now we check regime switching properly draws shocks
    mtest = AnSchorfheide(;testing = true)
    m = Model1002("ss10")
    m <= Setting(:saveroot, saveroot(mtest))
    m <= Setting(:dataroot, dataroot(mtest))
    m <= Setting(:regime_switching, true)
    m <= Setting(:use_parallel_workers, false)
    fp = dirname(@__FILE__)
    overrides = forecast_input_file_overrides(m)
    overrides[:full] = "$(fp)/../reference/mhsave_vint=181115.h5"
    θ = load_draws(m, :full)[1:3, :]
    m <= Setting(:forecast_jstep, 1)
    m <= Setting(:forecast_block_size, 3)

    rss = Vector{Dict{Int, Date}}(undef, 0)
    fss = Vector{Date}(undef, 0)
    ces = Vector{Date}(undef, 0)
    fcast_out = Dict()

    push!(rss, Dict{Int, Date}(1 => date_presample_start(m), 2 => Date(2020, 3, 31), 3 => Date(2020, 6, 30), 4 => Date(2020, 9, 30)))
    push!(fss, Date(2020, 6, 30))
    push!(ces, Date(2020, 6, 30))

    push!(rss, Dict{Int, Date}(1 => date_presample_start(m), 2 => Date(2020, 3, 31), 3 => Date(2020, 6, 30)))
    push!(fss, Date(2020, 6, 30))
    push!(ces, Date(2020, 6, 30))

    push!(rss, Dict{Int, Date}(1 => date_presample_start(m), 2 => Date(2020, 3, 31), 3 => Date(2020, 12, 31)))
    push!(fss, Date(2020, 6, 30))
    push!(ces, Date(2020, 12, 31))

    push!(rss, Dict{Int, Date}(1 => date_presample_start(m), 2 => Date(2020, 3, 31), 3 => Date(2021, 12, 31)))
    push!(fss, Date(2020, 6, 30))
    push!(ces, Date(2020, 12, 31))

    push!(rss, Dict{Int, Date}(1 => date_presample_start(m), 2 => Date(2020, 6, 30), 3 => Date(2021, 12, 31)))
    push!(fss, Date(2020, 6, 30))
    push!(ces, Date(2020, 12, 31))

    output_vars = [:forecastobs, :forecast4qobs, :forecastpseudo, :histpseudo]
    dfs = Dict()
    fp = dirname(@__FILE__)
    dfs[:full] = load(joinpath(fp, "../reference/regime_switch_data.jld2"), "regime_switch_df_full")
    dfs[:semi] = load(joinpath(fp, "../reference/regime_switch_data.jld2"), "regime_switch_df_semi")
    dfs[:none] = load(joinpath(fp, "../reference/regime_switch_data.jld2"), "regime_switch_df_none")

    if !generate_fulldist_forecast_data
        check_results = load(joinpath(fp, "../reference/regime_switching_fulldist_forecast_version=" * ver * ".jld2"), "fcast_out")
    end

    for (i, rs, fs, ce) in zip(1:length(rss), rss, fss, ces)
        # Set up regime switching
        m <= Setting(:regime_dates, rs)
        m <= Setting(:date_forecast_start, fs)
        m <= Setting(:date_conditional_end, ce)
        for i in 1:length(rs)
            adj = (i == 2) ? .9 : 1.
            ModelConstructors.set_regime_val!(m[:σ_g], i, adj * m[:σ_g].value)
            ModelConstructors.set_regime_val!(m[:σ_b], i, adj * m[:σ_b].value)
            ModelConstructors.set_regime_val!(m[:σ_μ], i, adj * m[:σ_μ].value)
            ModelConstructors.set_regime_val!(m[:σ_ztil], i, adj * m[:σ_ztil].value)
            ModelConstructors.set_regime_val!(m[:σ_λ_f], i, adj * m[:σ_λ_f].value)
            ModelConstructors.set_regime_val!(m[:σ_λ_w], i, adj * m[:σ_λ_w].value)
            ModelConstructors.set_regime_val!(m[:σ_r_m], i, adj * m[:σ_r_m].value)
            ModelConstructors.set_regime_val!(m[:σ_σ_ω], i, adj * m[:σ_σ_ω].value)
            ModelConstructors.set_regime_val!(m[:σ_π_star], i, adj * m[:σ_π_star].value)
            ModelConstructors.set_regime_val!(m[:σ_z_p], i, adj * m[:σ_z_p].value)
        end
        setup_regime_switching_inds!(m)

        # Construct parameter matrix
        θ_rs = zeros(size(θ, 1), ModelConstructors.n_parameters_regime_switching(m))
        n_pvec = DSGE.n_parameters(m)
        for i in 1:size(θ_rs, 1)
            θ_rs[i, 1:n_pvec] = map(x -> (abs(x) < 1e-14) ? 0. : x, θ[i, :])
            j = n_pvec
            for para in DSGE.get_parameters(m)
                if !isempty(para.regimes)
                    for (k, v) in para.regimes[:value]
                        if k > 1
                            j += 1
                            drawv = θ[i, m.keys[para.key]]
                            θ_rs[i, j] = (k  == 2) ? .9 * drawv : drawv
                            if abs(θ_rs[i, j]) < 1e-14
                                θ_rs[i, j] = 0.
                            end
                        end
                    end
                end
            end
        end

        fcast_out[i] = Dict()
        if date_conditional_end(m) == Date(2020, 12, 31) && i == 3
            lastrow = deepcopy(dfs[:full][end, :])
            push!(dfs[:full], vcat(Date(2020, 9, 30), lastrow[2:end]...))
            push!(dfs[:full], vcat(Date(2020, 12, 31), lastrow[2:end]...))
            lastrow = deepcopy(dfs[:semi][end, :])
            push!(dfs[:semi], vcat(Date(2020, 9, 30), lastrow[2:end]...))
            push!(dfs[:semi], vcat(Date(2020, 12, 31), lastrow[2:end]...))
        end
        for cond_type in [:full, :none, :semi]
            if i > 2 && cond_type == :none
                continue
            end
            forecast_one(m, :full, cond_type, output_vars, check_empty_columns = false,
                         params = θ_rs, df = dfs[cond_type], verbose = :none)
            output_files = get_forecast_output_files(m, :full, cond_type, output_vars)

            fcast_out[i][cond_type] = Dict()
            for var in keys(output_files)
                fcast_out[i][cond_type][var] = load(output_files[var], "arr")
            end

            if !generate_fulldist_forecast_data
                for var in keys(output_files)
                    @test maximum(abs.(check_results[i][cond_type][var] - fcast_out[i][cond_type][var])) < 5e-5
                end
            end
        end
    end

    if generate_fulldist_forecast_data
        fp = dirname(@__FILE__)
        JLD2.jldopen(joinpath(fp, "../reference/regime_switching_fulldist_forecast_version=" * ver * ".jld2"), true, true, true, IOStream) do file
            write(file, "fcast_out", fcast_out)
        end
    end
end

nothing


text:
#
# linalg.jl -
#
# A few linear algebra algorithms optimized conditionally.
#
module LinearAlgebraMethods

using ..MayOptimize
using ..MayOptimize: AVX, Standard

using LinearAlgebra
using LinearAlgebra: AbstractTriangular
import LinearAlgebra: dot, ldiv!, lmul!, cholesky, cholesky!

using Base: has_offset_axes
import Base: sum

abstract type AbstractAlgorithm end
abstract type CholeskyFactorization{opt<:OptimLevel} <: AbstractAlgorithm end
abstract type CholeskyLower{opt} <: CholeskyFactorization{opt} end
abstract type CholeskyUpper{opt} <: CholeskyFactorization{opt} end

# There are 2 variants to perform the Cholesky factorization
# Cholesky-Banachiewicz (which proceeds row-size) and Cholesky-Crout (which
# proceeds column-size) and they can be used to compute the L⋅L' or the R'⋅R
# factorizations.
struct CholeskyBanachiewiczLowerII{opt} <: CholeskyLower{opt} end
struct CholeskyBanachiewiczLowerI{ opt} <: CholeskyLower{opt} end
struct CholeskyBanachiewiczUpper{  opt} <: CholeskyUpper{opt} end
struct CholeskyCroutLower{         opt} <: CholeskyLower{opt} end
struct CholeskyCroutUpperII{       opt} <: CholeskyUpper{opt} end
struct CholeskyCroutUpperI{        opt} <: CholeskyUpper{opt} end

"""
     UntransposedLowerTriangular{T,S}

is the union of lower triangular matrix types defined in Julia standard library
`LinearAlgebra` and wrapped over regular matrices, that is
`LowerTriangular{T,S}` and `UnitLowerTriangular{T,S}` with `S` the type of an
untransposed matrix type whose elements have type `T`.  The matrix (of type
`S`) backing the storage of the entries is given by calling the `parent`
method.

For the parent, say `A`, of a lower-triangular matrix, only entries `A[i,j]`
for `i ≥ j` shall be accessed.

For the parent, say `A`, of a unit lower-triangular matrix, only entries
`A[i,j]` for `i > j` shall be accessed and `A[i,i] = 1` shall be assumed.

"""
const UntransposedLowerTriangular{T,S<:AbstractMatrix{T}} = Union{
    LowerTriangular{T,S},UnitLowerTriangular{T,S}}

"""
     UntransposedUpperTriangular{T,S}

is the union of upper triangular matrix types defined in standard Julia library
`LinearAlgebra` and wrapped over regular matrices, that is
`UpperTriangular{T,S}` and `UnitUpperTriangular{T,S}` with `S` the type of an
untransposed matrix type whose elements have type `T`.  The matrix (of type
`S`) backing the storage of the entries is given by calling the `parent`
method.

For the parent, say `A`, of an upper-triangular matrix, only entries `A[i,j]`
for `i ≤ j` shall be accessed.

For the parent, say `A`, of a unit upper-triangular matrix, only entries
`A[i,j]` for `i < j` shall be accessed and `A[i,i] = 1` shall be assumed.

"""
const UntransposedUpperTriangular{T,S<:AbstractMatrix{T}} = Union{
    UpperTriangular{T,S},UnitUpperTriangular{T,S}}

const TransposedLowerTriangular{T,S<:AbstractMatrix{T}} = Union{
    Adjoint{T,<:LowerTriangular{T,S}},
    Transpose{T,<:LowerTriangular{T,S}},
    Adjoint{T,<:UnitLowerTriangular{T,S}},
    Transpose{T,<:UnitLowerTriangular{T,S}},
    UpperTriangular{T,Adjoint{T,S}},
    UpperTriangular{T,Transpose{T,S}},
    UnitUpperTriangular{T,Adjoint{T,S}},
    UnitUpperTriangular{T,Transpose{T,S}}}

const TransposedUpperTriangular{T,S<:AbstractMatrix{T}} = Union{
    Adjoint{T,<:UpperTriangular{T,S}},
    Transpose{T,<:UpperTriangular{T,S}},
    Adjoint{T,<:UnitUpperTriangular{T,S}},
    Transpose{T,<:UnitUpperTriangular{T,S}},
    LowerTriangular{T,Adjoint{T,S}},
    LowerTriangular{T,Transpose{T,S}},
    UnitLowerTriangular{T,Adjoint{T,S}},
    UnitLowerTriangular{T,Transpose{T,S}}}

const TransposedTriangular{T,S<:AbstractMatrix{T}} = Union{
    TransposedLowerTriangular{T,S},
    TransposedUpperTriangular{T,S}}

const UntransposedTriangular{T,S<:AbstractMatrix{T}} = Union{
    UntransposedLowerTriangular{T,S},
    UntransposedUpperTriangular{T,S}}

default_optimization(::Type{<:AbstractAlgorithm}) = InBounds

for alg in (:CholeskyBanachiewiczLowerI,
            :CholeskyBanachiewiczLowerII,
            :CholeskyBanachiewiczUpper,
            :CholeskyCroutLower,
            :CholeskyCroutUpperI,
            :CholeskyCroutUpperII)
    @eval $alg(opt::Type{<:OptimLevel} = default_optimization($alg)) =
        $alg{opt}()
end

# Vectorization does not improve lower triangular Cholesky factorization, so
# make `InBounds` the default.
default_optimization(::Type{<:CholeskyLower}) = InBounds

# Vectorization does improve upper triangular Cholesky factorization, so make
# `Vectorize` the default.
default_optimization(::Type{<:CholeskyUpper}) = Vectorize

# Use the most efficient version by default.
const BestCholeskyFactorization = CholeskyBanachiewiczUpper
CholeskyFactorization() = BestCholeskyFactorization()
CholeskyFactorization(opt::Type{<:OptimLevel}) =
    BestCholeskyFactorization(opt)
CholeskyFactorization{opt}() where {opt<:OptimLevel} =
    BestCholeskyFactorization{opt}()

# Alias to the fastest Cholesky-Banachiewicz algorithm for the L'⋅L
# factorization.
const CholeskyBanachiewiczLower = CholeskyBanachiewiczLowerII

# Alias to the fastest Cholesky-Crout algorithm for the R⋅R' factorization.
const CholeskyCroutUpper = CholeskyCroutUpperII

const Floats = Union{AbstractFloat,Complex{<:AbstractFloat}}

#------------------------------------------------------------------------------
# Some basic linear algebra methods.

sum(::Type{<:Standard}, x::AbstractArray) = sum(x)

function sum(opt::Type{<:OptimLevel}, x::AbstractArray)
    s = zero(eltype(x))
    @maybe_vectorized opt for i in eachindex(x)
        s += x[i]
    end
    return s
end

dot(::Type{<:Standard}, x::AbstractVector, y::AbstractVector) = dot(x, y)

function dot(opt::Type{<:OptimLevel},
             x::AbstractArray{Tx,N},
             y::AbstractArray{Ty,N}) where {Tx,Ty,N}
    s = zero(promote_type(Tx, Ty))
    @maybe_vectorized opt for i in eachindex(x, y)
        s += conj(x[i])*y[i]
    end
    return s
end

#------------------------------------------------------------------------------
# Left division by a matrix.
#
# Call standard methods.
#
ldiv!(::Type{<:Standard}, A, B) = ldiv!(A, B)
ldiv!(::Type{<:Standard}, Y, A, B) = ldiv!(Y, A, B)
if VERSION < v"1.4.0-rc1"
    # Generic fallback. This assumes that B and Y have the same sizes.
    ldiv!(Y::AbstractArray, A::AbstractMatrix, B::AbstractArray) = begin
        Y === B || copyto!(Y, B)
        return ldiv!(A, Y)
    end
end

#
# Store A\b in y for A triangular.
#
function ldiv!(opt::Type{<:OptimLevel},
               y::AbstractVector{T},
               A::UntransposedTriangular{T},
               b::AbstractVector{T}) where {T<:Floats}
    y === b || copyto!(y, b)
    return ldiv!(opt, A, y)
end

#
# Store A\b in b for A upper triangular.
#
function ldiv!(opt::Type{<:OptimLevel},
               A::UntransposedUpperTriangular{T},
               b::AbstractVector{T}) where {T<:Floats}
    R = parent(A) # to avoid getindex overheads
    n = check_ldiv_args(R, b)
    @maybe_inbounds opt for j ∈ n:-1:1
        if b[j] != 0
            if !is_unit_triangular(A)
                b[j] /= R[j,j]
            end
            temp = b[j]
            @maybe_vectorized opt for i ∈ 1:j-1 # was j-1:-1:1 in BLAS
                b[i] -= temp*R[i,j]
            end
        end
    end
    return b
end
#
# Store A\b in b for A lower triangular.
#
function ldiv!(opt::Type{<:OptimLevel},
               A::UntransposedLowerTriangular{T},
               b::AbstractVector{T}) where {T<:Floats}
    L = parent(A) # to avoid getindex overheads
    n = check_ldiv_args(L, b)
    @maybe_inbounds opt for j ∈ 1:n
        if b[j] != 0
            if !is_unit_triangular(A)
                b[j] /= L[j,j]
            end
            temp = b[j]
            @maybe_vectorized opt for i ∈ j+1:n
                b[i] -= temp*L[i,j]
            end
        end
    end
    return b
end
#
# Store A'\b in b for A upper triangular.
#
function ldiv!(opt::Type{<:OptimLevel},
               A::TransposedUpperTriangular{T},
               b::AbstractVector{T}) where {T<:Floats}
    R = parent(parent(A)) # to avoid getindex overheads
    n = check_ldiv_args(R, b)
    f = elementwise_function(A)
    @maybe_inbounds opt for j ∈ 1:n
        temp = b[j]
        @maybe_vectorized opt for i ∈ 1:j-1
            temp -= f(R[i,j])*b[i]
        end
        if !is_unit_triangular(A)
            temp /= f(R[j,j])
        end
        b[j] = temp
    end
    return b
end

function ldiv!(opt::Type{<:OptimLevel},
               y::AbstractVector{T},
               A::TransposedUpperTriangular{T},
               b::AbstractVector{T}) where {T<:Floats}
    R = parent(parent(A)) # to avoid getindex overheads
    n = check_ldiv_args(y, R, b)
    f = elementwise_function(A)
    @maybe_inbounds opt for j ∈ 1:n
        temp = b[j]
        @maybe_vectorized opt for i ∈ 1:j-1
            temp -= f(R[i,j])*y[i]
        end
        if !is_unit_triangular(A)
            temp /= f(R[j,j])
        end
        y[j] = temp
    end
    return y
end
#
# Store A'\b in b for A lower triangular.
#
function ldiv!(opt::Type{<:OptimLevel},
               A::TransposedLowerTriangular{T},
               b::AbstractVector{T}) where {T<:Floats}
    L = parent(parent(A)) # to avoid getindex overheads
    n = check_ldiv_args(L, b)
    f = elementwise_function(A)
    @maybe_inbounds opt for j ∈ n:-1:1
        temp = b[j]
        @maybe_vectorized opt for i ∈ n:-1:j+1
            temp -= f(L[i,j])*b[i]
        end
        if !is_unit_triangular(A)
            temp /= f(L[j,j])
        end
        b[j] = temp
    end
    return b
end

function ldiv!(opt::Type{<:OptimLevel},
               y::AbstractVector{T},
               A::TransposedLowerTriangular{T},
               b::AbstractVector{T}) where {T<:Floats}
    L = parent(parent(A)) # to avoid getindex overheads
    n = check_ldiv_args(y, L, b)
    f = elementwise_function(A)
    @maybe_inbounds opt for j ∈ n:-1:1
        temp = b[j]
        @maybe_vectorized opt for i ∈ n:-1:j+1
            temp -= f(L[i,j])*y[i]
        end
        if !is_unit_triangular(A)
            temp /= f(L[j,j])
        end
        y[j] = temp
    end
    return y
end

function ldiv!(opt::Type{<:OptimLevel},
               A::Cholesky{T},
               b::AbstractVector{T}) where{T}
    uplo = getfield(A, :uplo)
    if uplo === 'U'
        R = getfield(A, :factors)
        ldiv!(opt, R, ldiv!(opt, R', b))
    elseif uplo === 'L'
        L = getfield(A, :factors)
        ldiv!(opt, L', ldiv!(opt, L, b))
    else
        throw_bad_uplo_field(A)
    end
    return b
end

function ldiv!(opt::Type{<:OptimLevel},
               y::AbstractVector{T},
               A::Cholesky{T},
               b::AbstractVector{T}) where{T}
    uplo = getfield(A, :uplo)
    if uplo === 'U'
        R = getfield(A, :factors)
        ldiv!(opt, R, ldiv!(opt, y, R', b))
    elseif uplo === 'L'
        L = getfield(A, :factors)
        ldiv!(opt, L', ldiv!(opt, y, L, b))
    else
        throw_bad_uplo_field(A)
    end
    return b
end

@noinline throw_bad_uplo_field(A::Cholesky) = error(
    "Unexpected field uplo='", getfield(A, :uplo),
    "' in Cholesky factorization")

#------------------------------------------------------------------------------

"""
    exec!(alg, dst, args...) -> dst

executes algorithm `alg` with arguments `args...` saving the result in
destination `dst`.

""" exec!

# Cholesky factorization can be done in-place.
exec!(alg::CholeskyFactorization, A::AbstractMatrix) = exec!(alg, A, A)

"""
    cholesky(alg, A) -> C

yields Cholesky factorization of matrix `A` using optimization level and/or
algorithm specified by `alg`.  For example, `alg` can be:

- `InBounds` or `Vectorize` to perform the fastest factorization without bounds
  checking or with loop vectorization respectively;

- `CholeskyBanachiewiczLower(Vectorize)` to use Cholesky-Banachiewicz algorithm
  with loop vectorization to compute the `L'⋅L` Cholesky factorization with `L`
  lower triangular;

- `CholeskyCroutUpper()` to use Cholesky-Crout algorithm with default
  optimization level to compute the `R⋅R'` Cholesky factorization with `R` upper
  triangular.

""" cholesky

"""
    cholesky!(alg, [buf,] A) -> C

yields Cholesky factorization of matrix `A` using optimization level and/or
algorithm specified by `alg`.  If argument `buf` is specified, the Cholesky
factorization is written in `buf` leaving `A` unchanged (unless `buf === A`);
otherwise, `A` is overwritten with the Cholesky factorization.

""" cholesky!

# Execute Cholesky method given optimization instance.
cholesky(opt::OptimLevel, A::AbstractMatrix) =
    cholesky(typeof(opt), A)
cholesky!(opt::OptimLevel, A::AbstractMatrix) =
    cholesky!(typeof(opt), A)
cholesky!(opt::OptimLevel, buf::AbstractMatrix, A::AbstractMatrix) =
    cholesky!(typeof(opt), buf, A)

# Call standard Cholesky method.
cholesky(::Type{<:Standard}, A::AbstractMatrix; kwds...) =
    cholesky(A, Val(false); kwds...)
cholesky!(::Type{<:Standard}, A::AbstractMatrix; kwds...) =
    cholesky!(A, Val(false); kwds...)
cholesky!(::Type{<:Standard}, buf::AbstractMatrix, A::AbstractMatrix; kwds...) =
    cholesky!(copyto!(buf, A), Val(false); kwds...)

# Call default Cholesky method given optimization type.
cholesky(opt::Type{<:OptimLevel}, A::AbstractMatrix) =
    cholesky(CholeskyFactorization(opt), A)
cholesky!(opt::Type{<:OptimLevel}, A::AbstractMatrix) =
    cholesky!(CholeskyFactorization(opt), A)
cholesky!(opt::Type{<:OptimLevel}, buf::AbstractMatrix, A::AbstractMatrix) =
    cholesky!(CholeskyFactorization{opt}, buf, A)

# Call specific Cholesky method given algorithm type.
cholesky(alg::Type{<:CholeskyFactorization}, A::AbstractMatrix) =
    cholesky(alg(), A)
cholesky!(alg::Type{<:CholeskyFactorization}, A::AbstractMatrix) =
    cholesky!(alg(), A)
function cholesky!(alg::Type{<:CholeskyFactorization}, buf::AbstractMatrix,
                   A::AbstractMatrix)
    return cholesky!(alg(), buf, A)
end

# Call specific Cholesky method given aklgorithm instance.
cholesky(alg::CholeskyFactorization, A::AbstractMatrix) =
    Cholesky(exec!(alg, similar(A), A), uplo_char(alg), 0)
cholesky!(alg::CholeskyFactorization, A::AbstractMatrix) =
    Cholesky(exec!(alg, A), uplo_char(alg), 0)
cholesky!(alg::CholeskyFactorization, buf::AbstractMatrix, A::AbstractMatrix) =
    Cholesky(exec!(alg, buf, A), uplo_char(alg), 0)

uplo_char(alg::CholeskyFactorization) = uplo_char(typeof(alg))
uplo_char(::Type{<:CholeskyLower}) = 'L'
uplo_char(::Type{<:CholeskyUpper}) = 'U'

#------------------------------------------------------------------------------
# Lower triangular Cholesky factorization.

# Cholesky–Banachiewicz algorithm (row-wise).
function exec!(alg::CholeskyBanachiewiczLowerI{opt},
               L::AbstractMatrix{T},
               A::AbstractMatrix{T}) where {T<:Floats,opt}
    n = check_chol_args(L, A)
    @maybe_inbounds opt for i ∈ 1:n
        @maybe_inbounds opt for j ∈ 1:i-1
            L[i,j] = (A[i,j] - _pdot(alg,L,i,j))/conj(L[j,j])
        end
        s = _pdot(alg,L,i)
        L[i,i] = sqrt(A[i,i] - s)
    end
    return L
end

# Improved version of Cholesky–Banachiewicz algorithm (row-wise).
function exec!(alg::CholeskyBanachiewiczLowerII{opt},
               L::AbstractMatrix{T},
               A::AbstractMatrix{T}) where {T<:Floats,opt}
    n = check_chol_args(L, A)
    @maybe_inbounds opt for i ∈ 1:n
        s = zero(real(T))
        @maybe_inbounds opt for j ∈ 1:i-1
            L[i,j] = (A[i,j] - _pdot(alg,L,i,j))/conj(L[j,j])
            s += abs2(L[i,j])
        end
        L[i,i] = sqrt(A[i,i] - s)
    end
    return L
end

# Cholesky–Crout algorithm (column-wise).
function exec!(alg::CholeskyCroutLower{opt},
               L::AbstractMatrix{T},
               A::AbstractMatrix{T}) where {T<:Floats,opt}
    n = check_chol_args(L, A)
    @maybe_inbounds opt for j ∈ 1:n
        s = _pdot(alg,L,j)
        L[j,j] = sqrt(A[j,j] - s)
        @maybe_inbounds opt for i ∈ j+1:n
            L[i,j] = (A[i,j] - _pdot(alg,L,i,j))/conj(L[j,j])
        end
    end
    return L
end

#------------------------------------------------------------------------------
# Upper triangular Cholesky factorization.

# Cholesky–Crout algorithm (column-wise).
function exec!(alg::CholeskyCroutUpperI{opt},
               R::AbstractMatrix{T},
               A::AbstractMatrix{T}) where {T<:Floats,opt}
    n = check_chol_args(R, A)
    @maybe_inbounds opt for j ∈ 1:n
        @maybe_inbounds opt for i ∈ 1:j-1
            R[i,j] = (A[i,j] - _pdot(alg,R,i,j))/conj(R[i,i])
        end
        s = _pdot(alg,R,j)
        R[j,j] = sqrt(A[j,j] - s)
    end
    return R
end

# Cholesky–Crout algorithm (column-wise) improved version.
function exec!(alg::CholeskyCroutUpperII{opt},
               R::AbstractMatrix{T},
               A::AbstractMatrix{T}) where {T<:Floats,opt}
    # Proceed column by column.
    n = check_chol_args(R, A)
    @maybe_inbounds opt for j ∈ 1:n
        s = zero(real(T))
        @maybe_inbounds opt for i ∈ 1:j-1
            R[i,j] = (A[i,j] - _pdot(alg,R,i,j))/conj(R[i,i])
            s += abs2(R[i,j])
        end
        R[j,j] = sqrt(A[j,j] - s)
    end
    return R
end

# Cholesky–Banachiewicz algorithm (row-wise).
function exec!(alg::CholeskyBanachiewiczUpper{opt},
               R::AbstractMatrix{T},
               A::AbstractMatrix{T}) where {T<:Floats,opt}
    n = check_chol_args(R, A)
    @maybe_inbounds opt for i ∈ 1:n
        s = _pdot(alg,R,i)
        R[i,i] = sqrt(A[i,i] - s)
        @maybe_inbounds opt for j ∈ i+1:n
            R[i,j] = (A[i,j] - _pdot(alg,R,i,j))/conj(R[i,i])
        end
    end
    return R
end

#------------------------------------------------------------------------------
# Partial row/column dot product for the Cholesky factorization.

"""
    _pdot(alg::CholeskyUpper, R, i, j) -> sum(conj(R[k,i])*R[k,j] for k ∈ 1:i-1)

yields the partial dot product of the `i`-th and `j`-th columns of matrix `R`
assuming it is a strict upper triangular matrix (as if entries are zero in the
diagonal and in the lower triangular part of `R`).

The condition `i < j` must hold, this is not checked.  If `i = j`, then call:

    _pdot(alg::CholeskyUpper, R, i) -> sum(abs2(R[k,i]) for k ∈ 1:i-1)

which is real even though `eltype(R)` may be complex.

---

    _pdot(alg::CholeskyLower, L, i, j) -> sum(L[i,k]*conj(L[j,k]) for k ∈ 1:j-1)

computes the partial dot product of the `i`-th and `j`-th rows of matrix `L`
assuming it is a strict lower triangular matrix (as if entries are zero in the
diagonal and in the upper triangular part of `L`).

The condition `i > j` must hold, this is not checked.  If `i = j`, then call:

    _pdot(alg::CholeskyLower, L, i) -> sum(abs2(L[i,k]) for k ∈ 1:i-1)

which is real even though `eltype(L)` may be complex.

---

The method `_pdot` is meant to be fast.  It is compiled according to the
optimization level embedded in `alg`.  The method shall be considered as
*unsafe* as, for maximum performances, the validity of the indices `i` and `j`
is not checked.

"""
@inline function _pdot(::CholeskyUpper{opt},
                       R::AbstractMatrix{T},
                       i::Int, j::Int) where {T,opt}
    s = zero(T)
    @maybe_vectorized opt for k ∈ 1:i-1
        s += conj(R[k,i])*R[k,j]
    end
    return s
end

@inline function _pdot(::CholeskyUpper{opt},
                       R::AbstractMatrix{T},
                       i::Int) where {T,opt}
    s = zero(real(T))
    @maybe_vectorized opt for k ∈ 1:i-1
        s += abs2(R[k,i])
    end
    return s
end

@inline function _pdot(::CholeskyLower{opt},
                       L::AbstractMatrix{T},
                       i::Int, j::Int) where {T,opt}
    s = zero(T)
    @maybe_vectorized opt for k in 1:j-1
        s += L[i,k]*conj(L[j,k])
    end
    return s
end

@inline function _pdot(::CholeskyLower{opt},
                       L::AbstractMatrix{T},
                       i::Int) where {T,opt}
    s = zero(real(T))
    @maybe_vectorized opt for k in 1:i-1
        s += abs2(L[i,k])
    end
    return s
end

#------------------------------------------------------------------------------

function check_chol_args(L::AbstractMatrix, A::AbstractMatrix)
    require_one_based_indexing(L, A)
    inds = axes(A)
    inds[1] == inds[2] || throw_dimension_mismatch(
        "expecting a square matrix")
    axes(L) == inds || throw_dimension_mismatch(
            "matrixes must have the same axes")
    return length(inds[1])
end

function check_ldiv_args(A::AbstractMatrix, b::AbstractVector)
    require_one_based_indexing(A, b)
    m, n = size(A)
    m == n || throw_dimension_mismatch(
        "A must be a square matrix")
    length(b) == n || throw_dimension_mismatch(
        "sizes of A and b are incompatible")
    return Int(n)
end

function check_ldiv_args(y::AbstractVector,
                         A::AbstractMatrix,
                         b::AbstractVector)
    require_one_based_indexing(y, A, b)
    m, n = size(A)
    m == n || throw_dimension_mismatch(
        "A must be a square matrix")
    length(y) == n || throw_dimension_mismatch(
        "sizes of A and y are incompatible")
    length(b) == n || throw_dimension_mismatch(
        "sizes of A and b are incompatible")
    return Int(n)
end

require_one_based_indexing(A...) =
    has_offset_axes(A...) && throw_argument_error(
        "arrays must have 1-based indexing")

throw_argument_error(msg::AbstractString) = throw(ArgumentError(msg))
@noinline throw_argument_error(args...) =
    throw_argument_error(string(args...))

throw_dimension_mismatch(msg::AbstractString) = throw(DimensionMismatch(msg))
@noinline throw_dimension_mismatch(args...) =
    throw_dimension_mismatch(string(args...))

"""
    is_unit_triangular(A)

yields whether the matrix `A` is a unit triangular matrix or not.  The result
is type-stable and is directly inferred from the type of `A`.  The argument can
also be a matrix type.

"""
is_unit_triangular(A::AbstractMatrix) = is_unit_triangular(typeof(A))
is_unit_triangular(::Type{<:Adjoint{T,S}}) where {T,S} = is_unit_triangular(S)
is_unit_triangular(::Type{<:Transpose{T,S}}) where {T,S} = is_unit_triangular(S)
is_unit_triangular(::Type{<:UnitLowerTriangular}) = true
is_unit_triangular(::Type{<:UnitUpperTriangular}) = true
is_unit_triangular(::Type{<:AbstractMatrix}) = false

"""
    elementwise_function(A)

yields `conj` if multiplying or dividing by the matrix `A` requires to take the
conjugate of the values of the object backing the storage of the entries of
`A`; otherwise yields `identity`.

"""
elementwise_function(A::AbstractMatrix) = elementwise_function(typeof(A))
elementwise_function(::Type{<:Adjoint}) = conj
elementwise_function(::Type{<:AbstractTriangular{<:Any,<:Adjoint}}) = conj
elementwise_function(::Type{<:AbstractMatrix}) = identity

end # module

ext: jl

max_stars_repo_name: emmt/ConditionallyOptimize


text:
module Test_pmapps

using Random
using LinearAlgebra
using MatrixPencils
using Polynomials
using Test

@testset "Polynomial Matrix Applications" begin

for fast in (false,true)

P = zeros(0,0,0);
@time @time val, kinfo  = pmeigvals(P, fast = fast)
@test val == [] && (kinfo.rki,kinfo.lki,kinfo.id,kinfo.nf,kinfo.nrank) == (Int64[], Int64[], Int64[], 0, 0)

P = zeros(0,0,3);
@time val, kinfo  = pmeigvals(P, fast = fast)
@test val == [] && (kinfo.rki,kinfo.lki,kinfo.id,kinfo.nf,kinfo.nrank) == (Int64[], Int64[], Int64[], 0, 0)

P = zeros(3,0,3);
@time val, kinfo  = pmeigvals(P, fast = fast)
@test val == [] && (kinfo.rki,kinfo.lki,kinfo.id,kinfo.nf,kinfo.nrank) == (Int64[], [0, 0, 0], Int64[], 0, 0)

P = zeros(0,3,3);
@time val, kinfo  = pmeigvals(P, fast = fast)
@test val == [] && (kinfo.rki,kinfo.lki,kinfo.id,kinfo.nf,kinfo.nrank) == ([0, 0, 0], Int64[], Int64[], 0, 0)

P = zeros(1,1,5)
@time val, kinfo  = pmeigvals(P, grade = pmdeg(P), fast = fast)
@test val == [] && (kinfo.rki,kinfo.lki,kinfo.id,kinfo.nf,kinfo.nrank) == (Int64[0], Int64[0], Int64[], 0, 0)

P = zeros(1,1,5)
@time val, kinfo  = pmeigvals(P, grade = pmdeg(P)+4, fast = fast)
@test val == [] && (kinfo.rki,kinfo.lki,kinfo.id,kinfo.nf,kinfo.nrank) == (Int64[0], Int64[0], Int64[], 0, 0)

P = zeros(1,1,5)
@time val, kinfo  = pmeigvals(P, CF1 = false, grade = pmdeg(P)+3, fast = fast)
@test val == [] && (kinfo.rki,kinfo.lki,kinfo.id,kinfo.nf,kinfo.nrank) == (Int64[0], Int64[0], Int64[], 0, 0)

P = zeros(1,1,5)
P[1,1,1] = 1.

@time val, kinfo  = pmeigvals(P, fast = fast)
@test val == Float64[] && (kinfo.rki,kinfo.lki,kinfo.id,kinfo.nf,kinfo.nrank) == (Int64[], Int64[], Int64[], 0, 1)

@time val, kinfo  = pmeigvals(P, grade = 1, fast = fast)
@test val == [Inf] && (kinfo.rki,kinfo.lki,kinfo.id,kinfo.nf,kinfo.nrank) == (Int64[], Int64[], Int64[1], 0, 1)

@time val, kinfo  = pmeigvals(P, grade = pmdeg(P)+4, fast = fast)
@test val == [Inf, Inf, Inf, Inf] && (kinfo.rki,kinfo.lki,kinfo.id,kinfo.nf,kinfo.nrank) == (Int64[], Int64[], Int64[4], 0, 1)

@time val, kinfo  = pmeigvals(P, CF1 = false, grade = pmdeg(P)+3, fast = fast)
@test val == [Inf, Inf, Inf] && (kinfo.rki,kinfo.lki,kinfo.id,kinfo.nf,kinfo.nrank) == (Int64[], Int64[], Int64[3], 0, 1)

P = zeros(1,1,5)
P[1,1,3] = 1.
P[1,1,2] = -3.
P[1,1,1] = 2.
@time val, kinfo  = pmeigvals(P, grade = pmdeg(P), fast = fast)
@test (val ≈ [2, 1] || val ≈ [1, 2]) && (kinfo.rki,kinfo.lki,kinfo.id,kinfo.nf,kinfo.nrank) == (Int64[], Int64[], Int64[], 2, 1) &&
kinfo.nf+sum(kinfo.id)+sum(kinfo.rki)+sum(kinfo.lki) == kinfo.nrank*pmdeg(P)

@time val, kinfo  = pmeigvals(P, grade = 4, fast = fast)
@test sort(val) ≈ [1, 2, Inf, Inf] && (kinfo.rki,kinfo.lki,kinfo.id,kinfo.nf,kinfo.nrank) == (Int64[], Int64[], Int64[2], 2, 1) &&
      kinfo.nf+sum(kinfo.id)+sum(kinfo.rki)+sum(kinfo.lki) == kinfo.nrank*4

@time val, kinfo  = pmeigvals(P, CF1 = false, grade = 3, fast = fast)
@test sort(val) ≈ [1, 2, Inf] && (kinfo.rki,kinfo.lki,kinfo.id,kinfo.nf,kinfo.nrank) == (Int64[], Int64[], Int64[1], 2, 1) &&
kinfo.nf+sum(kinfo.id)+sum(kinfo.rki)+sum(kinfo.lki) == kinfo.nrank*3

# P = 0

@time @time info, iz, ip = pmkstruct(0, fast = fast)
@test (info.rki, info.lki,info.id, info.nf, info.nrank) == ([0], [0], [], 0, 0) && iz == [] && ip == []

@time @time info, iz, ip = pmkstruct(0,grade = 3, fast = fast)
@test (info.rki, info.lki,info.id, info.nf, info.nrank) == ([0], [0], [], 0, 0) && iz == [] && ip == []

@time @time val, info = pmeigvals(0, fast = fast)
@test val ≈ Float64[] && (info.rki, info.lki,info.id, info.nf,info.nrank) == ([0], [0], [], 0, 0)

@time val, info = pmeigvals(0,grade=1, fast = fast)
@test val ≈ Float64[] && (info.rki, info.lki,info.id, info.nf,info.nrank) == ([0], [0], [], 0, 0)

@time val, info = pmeigvals(0,grade=3, fast = fast)
@test val ≈ Float64[] && (info.rki, info.lki,info.id, info.nf,info.nrank) == ([0], [0], [], 0, 0)

@time val, iz, info = pmzeros(0, fast = fast)
@test val ≈ Float64[] && iz == [] && (info.rki, info.lki,info.id, info.nf,info.nrank) == ([0], [0], [], 0, 0)

@time val, ip, id = pmpoles(0, fast = fast)
@test val ≈ Float64[] && ip == [] && id == []

@test pmrank(0) == 0 && pmrank(0,fastrank=false) == 0

@test !ispmregular(0) && !ispmregular(0,fastrank=false)

@test !ispmunimodular(0)

# P = 1

@time @time info, iz, ip = pmkstruct(1, fast = fast)
@test (info.rki, info.lki,info.id, info.nf, info.nrank) == ([], [], [], 0, 1) && iz == [] && ip == []

@time @time info, iz, ip = pmkstruct(1,grade=3, fast = fast)
@test (info.rki, info.lki,info.id, info.nf, info.nrank) == ([], [], [3], 0, 1) && iz == [] && ip == []

@time val, info = pmeigvals(1, fast = fast)
@test val ≈ Float64[] && (info.rki, info.lki,info.id, info.nf,info.nrank) == ([], [], [], 0, 1)

@time val, info = pmeigvals(1,grade=1, fast = fast)
@test val ≈ [Inf] && (info.rki, info.lki,info.id, info.nf,info.nrank) == ([], [], [1], 0, 1)

@time val, info = pmeigvals(1,grade=3, fast = fast)
@test val ≈ [Inf, Inf, Inf] && (info.rki, info.lki,info.id, info.nf,info.nrank) == ([], [], [3], 0, 1)

@time val, iz, info = pmzeros(1, fast = fast)
@test val ≈ Float64[] && iz == [] && (info.rki, info.lki,info.id, info.nf,info.nrank) == ([], [], [], 0, 1)

@time val, ip, id = pmpoles(1, fast = fast)
@test val ≈ Float64[] && ip == [] && id == []

@test pmrank(1) == 1 && pmrank(1,fastrank=false) == 1

@test ispmregular(1) && ispmregular(1,fastrank=false)

@test ispmunimodular(1)

# P = λ

λ = Polynomial([0,1],:λ)

@time @time info, iz, ip = pmkstruct(λ, fast = fast)
@test (info.rki, info.lki,info.id, info.nf, info.nrank) == ([], [], [], 1, 1) && iz == [] && ip == [1]

@time @time info, iz, ip = pmkstruct(λ,grade=3, fast = fast)
@test (info.rki, info.lki,info.id, info.nf, info.nrank) == ([], [], [2], 1, 1) && iz == [] && ip == [1]

@time val, info = pmeigvals(λ, fast = fast)
@test val ≈ [0.] && (info.rki, info.lki,info.id, info.nf,info.nrank) == ([], [], [], 1, 1)

@time val, info = pmeigvals([λ 1;1 0], fast = fast)
@test val ≈ [Inf, Inf] && (info.rki, info.lki,info.id, info.nf,info.nrank) == ([], [], [2], 0, 2)

@time val, info = pmeigvals(λ,grade=3, fast = fast)
@test val ≈ [0., Inf, Inf] && (info.rki, info.lki,info.id, info.nf,info.nrank) == ([], [], [2], 1, 1)

@time val, info = pmeigvals([λ 1;1 0],grade=3, fast = fast)
@test val ≈ [Inf, Inf, Inf, Inf, Inf, Inf] && (info.rki, info.lki,info.id, info.nf,info.nrank) == ([], [], [2, 4], 0, 2)

@time val, iz, info = pmzeros(λ, fast = fast)
@test val ≈ [0.] && iz == [] && (info.rki, info.lki,info.id, info.nf,info.nrank) == ([], [], [], 1, 1)

val = pmroots(λ, fast = fast)
@test val ≈ [0.]

@time @time val, ip, id = pmpoles(λ, fast = fast)
@test val == [Inf] && ip == [1] && id == [2]

@test pmrank(λ) == 1 && pmrank(λ,fastrank=false) == 1

@test ispmregular(λ) && ispmregular(λ,fastrank=false)

@test !ispmunimodular(λ)

# P = [λ 1]
λ = Polynomial([0,1],:λ)
P = [λ one(λ)]

@time @time info, iz, ip = pmkstruct(P, fast = fast)
@test (info.rki, info.lki,info.id, info.nf, info.nrank) == ([1], [], [], 0, 1) && iz == [] && ip == [1]

@time @time info, iz, ip = pmkstruct(P,CF1=true, fast = fast)
@test (info.rki, info.lki,info.id, info.nf, info.nrank) == ([1], [], [], 0, 1) && iz == [] && ip == [1]

@time @time info, iz, ip = pmkstruct(P,grade=3, fast = fast)
@test (info.rki, info.lki,info.id, info.nf, info.nrank) == ([1], [], [2], 0, 1) && iz == [] && ip == [1]

@time @time info, iz, ip = pmkstruct(P,grade=3,CF1=true, fast = fast)
@test (info.rki, info.lki,info.id, info.nf, info.nrank) == ([1], [], [2], 0, 1) && iz == [] && ip == [1]

@time val, info = pmeigvals(P, fast = fast)
@test val == Float64[] && (info.rki, info.lki,info.id, info.nf,info.nrank) == ([1], [], [], 0, 1)

@time val, info = pmeigvals(P,grade=3, fast = fast)
@test val == [Inf, Inf] && (info.rki, info.lki,info.id, info.nf,info.nrank) == ([1], [], [2], 0, 1)

@time val, iz, info = pmzeros(P, fast = fast)
@test val == Float64[] && iz == [] && (info.rki, info.lki,info.id, info.nf,info.nrank) == ([1], [], [], 0, 1)

@time val, ip, id = pmpoles(P, fast = fast)
@test val == [Inf] && ip == [1] && id == [1, 2]

@test pmrank(P) == 1 && pmrank(P,fastrank=false) == 1

@test !ispmregular(P) && !ispmregular(P,fastrank=false)

@test !ispmunimodular(P)

# Example 3: P = [λ^2 λ; λ 1] DeTeran, Dopico, Mackey, ELA 2009
λ = Polynomial([0,1],:λ)
P = [λ^2 λ; λ 1]

@time info, iz, ip = pmkstruct(P, fast = fast)
@test (info.rki, info.lki,info.id, info.nf, info.nrank) == ([1], [1], [], 0, 1) && iz == [] && ip == [2]

@time info, iz, ip = pmkstruct(P,CF1=true, fast = fast)
@test (info.rki, info.lki,info.id, info.nf, info.nrank) == ([1], [1], [], 0, 1) && iz == [] && ip == [2]

@time info, iz, ip = pmkstruct(P,grade=3, fast = fast)
@test (info.rki, info.lki,info.id, info.nf, info.nrank) == ([1], [1], [1], 0, 1) && iz == [] && ip == [2]

@time info, iz, ip = pmkstruct(P,grade=3,CF1=true, fast = fast)
@test (info.rki, info.lki,info.id, info.nf, info.nrank) == ([1], [1], [1], 0, 1) && iz == [] && ip == [2]

@time val, info = pmeigvals(P, fast = fast)
@test val == Float64[] && (info.rki, info.lki,info.id, info.nf,info.nrank) == ([1], [1], [], 0, 1)

@time val, info = pmeigvals(P,grade=3, fast = fast)
@test val == [Inf] && (info.rki, info.lki,info.id, info.nf,info.nrank) == ([1], [1], [1], 0, 1)

@time val, iz, info = pmzeros(P, fast = fast)
@test val == Float64[] && iz == [] && (info.rki, info.lki,info.id, info.nf,info.nrank) == ([1], [1], [], 0, 1)

@time val, ip, id = pmpoles(P, fast = fast)
@test val == [Inf, Inf] && ip == [2] && id == [2, 2, 4]

@test pmrank(P) == 1 && pmrank(P,fastrank=false) == 1

@test !ispmregular(P) && !ispmregular(P,fastrank=false)

@test !ispmunimodular(P)

# Example: DeTeran, Dopico, Mackey, ELA 2009
P = zeros(2,2,3);
P[:,:,1] = [0 0; 0 1.];
P[:,:,2] = [0 1.; 1. 0];
P[:,:,3] = [1. 0; 0 0];
@time val, info = pmeigvals(P)
@test val == Float64[] && (info.rki, info.lki,info.id, info.nf,info.nrank) == ([1], [1], [], 0, 1)

# Example 4: Van Dooren, Dewilde, LAA, 1983
P = zeros(Int,3,3,3)
P[:,:,1] = [1 2 -2; 0 -1 -2; 0 0 0]
P[:,:,2] = [1 3 0; 1 4 2; 0 -1 -2]
P[:,:,3] = [1 4 2; 0 0 0; 1 4 2]

@time val, info = pmeigvals(P, fast = fast)
@test val ≈ [1, Inf, Inf] && (info.rki, info.lki,info.id, info.nf,info.nrank) == ([0], [1], [2], 1, 2)

@time val, info = pmeigvals(P,grade = 3, fast = fast)
@test val ≈ [1, Inf, Inf, Inf, Inf] && (info.rki, info.lki,info.id, info.nf,info.nrank) == ([0], [1], [1, 3], 1, 2)

@time val, iz, info = pmzeros(P, fast = fast)
@test val ≈ [1] && iz == [] && (info.rki, info.lki,info.id, info.nf,info.nrank) == ([0], [1], [2], 1, 2)

@time val, ip, id = pmpoles(P, fast = fast)
@test val ≈ [Inf, Inf] && ip == [2] && id == [2, 2, 2, 2, 4]

val, iz, info  = pmzeros1(P, fast = fast)
@test val ≈ [1] && iz == [] && (info.rki, info.lki,info.id, info.nf,info.nrank) == ([0], [1], [1], 1, 2)

@time val, ip, id = pmpoles1(P, fast = fast)
@test val ≈ [Inf, Inf] && ip == [2] && id == [1, 1, 1, 1, 3]

@time val, iz, info = pmzeros2(P, fast = fast)
@test val ≈ [1] && iz == [] && (info.rki, info.lki,info.id, info.nf,info.nrank) == ([0], [1], [1, 1, 1, 1], 1, 2)

@time val, ip, id = pmpoles2(P, fast = fast)
@test val ≈ [Inf, Inf] && ip == [2] && id == [1, 3]

@test !ispmregular(P) && !ispmregular(P,fastrank=false)

@test !ispmunimodular(P)

λ = Polynomial([0,1],:λ);
P = [λ^2 + λ + 1 4λ^2 + 3λ + 2 2λ^2 - 2;
λ 4λ - 1 2λ - 2;
λ^2 4λ^2 - λ 2λ^2 - 2λ]

@time info, iz, ip = pmkstruct(P, fast = fast)
@test (info.rki, info.lki,info.id, info.nf, info.nrank) == ([0], [1], [2], 1, 2)  && iz == [] && ip == [2]

@time info, iz, ip = pmkstruct(P,CF1=true, fast = fast)
@test (info.rki, info.lki,info.id, info.nf, info.nrank) == ([0], [1], [2], 1, 2)  && iz == [] && ip == [2]

@time info, iz, ip = pmkstruct(P,grade=3, fast = fast)
@test (info.rki, info.lki,info.id, info.nf, info.nrank) == ([0], [1], [1, 3], 1, 2)  && iz == [] && ip == [2]

@time info, iz, ip = pmkstruct(P,grade=3,CF1=true, fast = fast)
@test (info.rki, info.lki,info.id, info.nf, info.nrank) == ([0], [1], [1, 3], 1, 2)  && iz == [] && ip == [2]

@time val, info = pmeigvals(P, fast = fast)
@test val ≈ [1, Inf, Inf] && (info.rki, info.lki,info.id, info.nf,info.nrank) == ([0], [1], [2], 1, 2)

@time val, info = pmeigvals(P,grade = 3, fast = fast)
@test val ≈ [1, Inf, Inf, Inf, Inf] && (info.rki, info.lki,info.id, info.nf,info.nrank) == ([0], [1], [1, 3], 1, 2)

@time val, info = pmeigvals(P,CF1=false, grade = 3, fast = fast)
@test val ≈ [1, Inf, Inf, Inf, Inf] && (info.rki, info.lki,info.id, info.nf,info.nrank) == ([0], [1], [1, 3], 1, 2)

@time val, iz, info = pmzeros(P, fast = fast)
@test val ≈ [1] && iz == [] && (info.rki, info.lki,info.id, info.nf,info.nrank) == ([0], [1], [2], 1, 2)

@time val, iz, id = pmpoles(P, fast = fast)
@test val ≈ [Inf, Inf] && iz == [2] && id == [2, 2, 2, 2, 4]

@time val, iz, id = pmpoles(P,CF1=false, fast = fast)
@test val ≈ [Inf, Inf] && iz == [2] && id == [2, 2, 2, 2, 4]

@time val, iz, info = pmzeros1(P, fast = fast)
@test val ≈ [1] && iz == [] && (info.rki, info.lki,info.id, info.nf,info.nrank) == ([0], [1], [1], 1, 2)

@time val, ip, id = pmpoles1(P, fast = fast)
@test val ≈ [Inf, Inf] && ip == [2] && id == [1, 1, 1, 1, 3]

@time val, iz, info = pmzeros2(P, fast = fast)
@test val ≈ [1] && iz == [] && (info.rki, info.lki,info.id, info.nf,info.nrank) == ([0], [1], [1, 1, 1, 1], 1, 2)

@time val, ip, id = pmpoles2(P, fast = fast)
@test val ≈ [Inf, Inf] && ip == [2] && id == [1, 3]

@test pmrank(P) == 2 && pmrank(P,fastrank=false) == 2

@test pmrank(rand(2,3,5)) == 2 && pmrank(rand(2,3,5),fastrank=false) == 2

@test pmrank(rand(3,2,5)) == 2 && pmrank(rand(3,2,5),fastrank=false) == 2

@test !ispmregular(P) && !ispmregular(P,fastrank=false)

@test !ispmunimodular(P)

s = Polynomial([0, 1],:x)
P = [ s^2+s+1  4*s^2+3*s+2 2*s^2-2
      s  4*s-1  2*s-2
      s^2 4*s^2-s  2*s^2-2*s ]
U = [-s  0   s^2 + 2*s + 1
s^2 - s + 1  0 -s^3 - s^2 - 1
0 -1 s];
V = [3*s + 1 3*s^3 + s^2 - 3*s - 4  6
-s -s^3 + s + 1  -2
0  0  1];
@test ispmunimodular(U)
@test ispmunimodular(V)

# unimodular matrix example
P = zeros(2,2,2)
P[:,:,1] = [1 0;0 1]; P[:,:,2] = [0. 0; 1 0];
@test ispmunimodular(P)
@time val, kinfo  = pmeigvals(P, CF1 = true, fast = fast)
@test val == [Inf, Inf] && (kinfo.rki,kinfo.lki,kinfo.id,kinfo.nf,kinfo.nrank) == (Int64[], Int64[], Int64[2], 0, 2)
@time val, kinfo  = pmeigvals(P, CF1 = false, fast = fast)
@test val == [Inf, Inf] && (kinfo.rki,kinfo.lki,kinfo.id,kinfo.nf,kinfo.nrank) == (Int64[], Int64[], Int64[2], 0, 2)
val, iz, kinfo  = pmzeros(P, CF1 = true, fast = fast)
@test val == [Inf] && iz == [1] && (kinfo.rki,kinfo.lki,kinfo.id,kinfo.nf,kinfo.nrank) == (Int64[], Int64[], Int64[2], 0, 2)
val, iz, kinfo  = pmzeros(P, CF1 = false, fast = fast)
@test val == [Inf] && iz == [1] && (kinfo.rki,kinfo.lki,kinfo.id,kinfo.nf,kinfo.nrank) == (Int64[], Int64[], Int64[2], 0, 2)
val, iz, kinfo  = pmzeros2(P, fast = fast)
@test val == [Inf] && iz == [1] && (kinfo.rki,kinfo.lki,kinfo.id,kinfo.nf,kinfo.nrank) == (Int64[], Int64[], Int64[1, 1, 2], 0, 2)

N = zeros(1,1); M = ones(1,1);
P = zeros(1,1,3)
P[:,:,1] = M; P[:,:,2] = -N;
val, iz, kinfo  = pmzeros2(P, fast = fast)
@test val == [] && iz == [] && (kinfo.rki,kinfo.lki,kinfo.id,kinfo.nf,kinfo.nrank) == (Int64[], Int64[], Int64[1], 0, 1)

M = [1 0;0 1]; N = [0 1; 0 0];
P = zeros(2,2,3)
P[:,:,1] = M; P[:,:,2] = -N;
val, kinfo = pmeigvals(P, fast = fast)
@test val == [Inf;Inf] && (kinfo.rki,kinfo.lki,kinfo.id,kinfo.nf,kinfo.nrank) == (Int64[], Int64[], [2], 0, 2)

M = [  22  34  31   31  17
        45  45  42   19  29
        39  47  49   26  34
        27  31  26   21  15
        38  44  44   24  30 ];

N = [   13  26  25  17  24
        31  46  40  26  37
        26  40  19  25  25
        16  25  27  14  23
        24  35  18  21  22  ];

P = zeros(5,5,3)
P[:,:,1] = M; P[:,:,2] = -N;
val, kinfo = pmeigvals(P, fast = fast)
@test (isapprox(val,[2,0,0,Inf],atol = 1.e-7) || isapprox(val,[0,0,2,Inf],atol = 1.e-7)) &&
      (kinfo.rki,kinfo.lki,kinfo.id,kinfo.nf,kinfo.nrank) == ([0], [0], [1], 3, 4)

# Test Suite 2 (Example 2.2.1, Beelen)
M = zeros(14,16);
i1 = [1 2 3 6 7 8 9 10 11 12 13 13 14];
j1 = [4 6 7 8 9 10 11 12 13 14 15 16 16];
aval = [1 1 1 1 1 1 1 1 1 2 3 1 3];
for k = 1:length(i1)
    M[i1[k],j1[k]] = aval[k]
end
N = zeros(14,16);
i1 = [1 2 3 5 6 7 10 12 13 14]; j1 = [3 5 6 8 9 10 13 14 15 16];
for k = 1:length(i1)
    N[i1[k],j1[k]] = 1
end

Q = qr(rand(14,14)).Q; Z = qr(rand(16,16)).Q;
M2 = Q*M*Z; N2 = Q*N*Z;

M = copy(M2); N = copy(N2);
P = zeros(size(M)...,3)
P[:,:,1] = M; P[:,:,2] = -N;

val, kinfo = pmeigvals(P, fast = fast, atol = 1.e-7)
@test length(val) == 6 &&
      length(filter(y-> y == true,isinf.(val))) == 3 &&
      length(filter(y-> y == true,isfinite.(val))) == 3 &&
      (kinfo.rki,kinfo.lki,kinfo.id,kinfo.nf,kinfo.nrank) == ([0, 0, 1, 2], [0, 3], [1, 2], 3, 12)

for Ty in (Float64, Complex{Float64})

P = rand(Ty,2,3,5);
P1 = rand(Ty,3,2,5);
abstol = sqrt(eps(one(real(Ty))))

@time info, iz, ip = pmkstruct(P,CF1=false, fast = fast, atol = abstol)
@test (info.rki, info.lki,info.id, info.nf, info.nrank) == ([8], [], [], 0, 2) && iz == [] && ip ==[4, 4]

@time info, iz, ip = pmkstruct(P,CF1=true, fast = fast, atol = abstol)
@test (info.rki, info.lki,info.id, info.nf, info.nrank) == ([8], [], [], 0, 2) && iz == [] && ip ==[4, 4]

@time info, iz, ip = pmkstruct(P,grade=5, CF1=false, fast = fast, atol = abstol)
@test (info.rki, info.lki,info.id, info.nf, info.nrank) == ([8], [], [1, 1], 0, 2) && iz == [] && ip ==[4, 4]

@time info, iz, ip = pmkstruct(P,grade=5,CF1=true, fast = fast, atol = abstol)
@test (info.rki, info.lki,info.id, info.nf, info.nrank) == ([8], [], [1, 1], 0, 2) && iz == [] && ip ==[4, 4]

@time val, info = pmeigvals(P, fast = fast, atol = abstol)
@test val == Ty[] && (info.rki, info.lki,info.id, info.nf,info.nrank) == ([8], [], [], 0, 2)

@time val, info = pmeigvals(P,CF1=true,grade = 5, fast = fast, atol = abstol)
@test val == [Inf, Inf] && (info.rki, info.lki,info.id, info.nf,info.nrank) == ([8], [], [1, 1], 0, 2)

@time val, info = pmeigvals(P,CF1=false, grade = 5, fast = fast, atol = abstol)
@test val == [Inf, Inf] && (info.rki, info.lki,info.id, info.nf,info.nrank) == ([8], [], [1, 1], 0, 2)

@time val, iz, info = pmzeros(P, fast = fast, atol = abstol)
@test val == Ty[] && iz == [] && (info.rki, info.lki,info.id, info.nf,info.nrank) == ([8], [], [], 0, 2)

@time val, ip, id = pmpoles(P,CF1=false, fast = fast, atol = abstol)
@test val == [Inf, Inf, Inf, Inf, Inf, Inf, Inf, Inf] && ip == [4, 4] && id == [4, 8, 8]

@time val, ip, id = pmpoles(P,CF1=true, fast = fast, atol = abstol)
@test val == [Inf, Inf, Inf, Inf, Inf, Inf, Inf, Inf] && ip == [4, 4] && id == [4, 8, 8]

@time val, iz, info = pmzeros1(P, fast = fast, atol = abstol)
@test val == Ty[] && iz == [] && (info.rki, info.lki,info.id, info.nf,info.nrank) == ([8], [], [], 0, 2)

@time val, iz, info = pmzeros1(P1, fast = fast, atol = abstol)
@test val == Ty[] && iz == [] && (info.rki, info.lki,info.id, info.nf,info.nrank) == ([], [8], [], 0, 2)

@time val, ip, id = pmpoles1(P, fast = fast, atol = abstol)
@test val == [Inf, Inf, Inf, Inf, Inf, Inf, Inf, Inf] && ip == [4, 4] && id == [1, 5, 5]

@time val, ip, id = pmpoles1(P1, fast = fast, atol = abstol)
@test val == [Inf, Inf, Inf, Inf, Inf, Inf, Inf, Inf] && ip == [4, 4] && id == [1, 5, 5]

@time val, iz, info = pmzeros2(P, fast = fast, atol = abstol)
@test val == Ty[] && iz == [] && (info.rki, info.lki,info.id, info.nf,info.nrank) == ([8], [], [1, 1, 1, 1], 0, 2)

@time val, ip, id = pmpoles2(P, fast = fast, atol = abstol)
@test val == [Inf, Inf, Inf, Inf, Inf, Inf, Inf, Inf] && ip == [4, 4] && id == [5, 5]

@test pmrank(P) == 2 && pmrank(P,fastrank=false) == 2

@test !ispmregular(P) && !ispmregular(P,fastrank=false)

@test !ispmunimodular(P)

@test ispmregular(rand(Ty,2,2,5)) && ispmregular(rand(Ty,2,2,5),fastrank=false)

@test !ispmunimodular(rand(Ty,2,2,5))

@test ispmunimodular(rand(Ty,2,2,1))

end
end
#end

end  # applications testset

end  # module