Jutho

# gkl.jl
"""
    mutable struct GKLFactorization{TU,TV,S<:Real}

Structure to store a Golub-Kahan-Lanczos (GKL) bidiagonal factorization of a linear map `A`
of the form

```julia
A * V = U * B + r * b'
A' * U = V * B'

Jutho

; Function rec_diff
  %2 = getelementptr [20 x i64], [20 x i64] addrspace(11)* %1, i64 0, i64 1
  %3 = getelementptr [20 x i64], [20 x i64] addrspace(11)* %1, i64 0, i64 0
  %4 = load i64, i64 addrspace(11)* %3, align 8
  %5 = getelementptr [20 x i64], [20 x i64] addrspace(11)* %1, i64 0, i64 5
  %6 = getelementptr [20 x i64], [20 x i64] addrspace(11)* %1, i64 0, i64 9
  %7 = getelementptr [20 x i64], [20 x i64] addrspace(11)* %1, i64 0, i64 13
  %8 = getelementptr [20 x i64], [20 x i64] addrspace(11)* %1, i64 0, i64 17
  %9 = getelementptr [20 x i64], [20 x i64] addrspace(11)* %1, i64 0, i64 18
  %10 = bitcast i64 addrspace(11)* %2 to <4 x i64> addrspace(11)*

Jutho

import Base.LinAlg: SingularException, checksquare

# QR ALGORITHM:
# Compute Schur form of a Hessenberg matrix

# Givens transform the upper right corner of a matrix.
function transform!(H::AbstractMatrix, U, G::Givens, imax::Int = size(H,1), jmin::Int = 1)
    lmul!(H, G, jmin:size(H,2))
    rmulc!(H, G, 1:imax)
    rmulc!(U, G)

Jutho

src/transpose.cpp:383:39: error: assigning to 'float' from incompatible type
      'complex<float>'
               B[i + j * ldb] = alpha * std::conj(A[i * lda + j]);
                                ~~~~~~^~~~~~~~~~~~~~~~~~~~~~~~~~~
src/transpose.cpp:676:13: note: in instantiation of function template
      specialization 'hptt::macro_kernel_scalar<1, float, true>' requested here
            macro_kernel_scalar<betaIsZero,floatType, conjA>(&A[i*lda], ...
            ^
src/transpose.cpp:842:16: note: in instantiation of function template
      specialization 'hptt::transpose_int<32, 32, 1, float, false, true>'

Jutho

import Base: start, next, done, size, length, eltype, iteratorsize, iteratoreltype, first, last
using Base: tail, IteratorEltype, HasEltype, EltypeUnknown, IteratorSize, HasShape, HasLength, IsInfinite, SizeUnknown
immutable ProductIterator{T<:Tuple}
    iterators::T
end

ProductIterator(iterators...) = ProductIterator(iterators)

iteratorsize(P::ProductIterator) = _itsize(P.iterators)
@inline _itsize(::Tuple{}) = HasShape()

Jutho

Preliminary notes:

• Abstract is left out of the type names for conciseness
• We need a Covector (or alternative name) type for theoretical reasons, i.e. to make dispatch able to distinguish between all of the following operations and to have methods with a different return type for these different operations.
• We need a LazyTranspose type (and Covector) for computational efficiency, i.e. to link with BLAS, to replace the current Ac_mul_Bt etc family.
• Alternatives with ? ... ? are not very likely, just to point out the possibility.

Given vectors v and w, linear maps with matrix representation A and B, and a covector z

abstract operation	julia ascii	julia method	possible alternative or unicode equivalent

Jutho

function analyzeconvert()
    methodlist = methods(convert,Tuple{Type,Number})

    dict = Dict{DataType,Vector{Any}}()
    for m in methodlist
        to = m.sig.parameters[1].parameters[1]
        from = m.sig.parameters[2]

        TT = isa(to,Union) ? to.types : Base.svec(to)
        TF = isa(from,Union) ? from.types : Base.svec(from)

Jutho

### Multidimensional iterators
module IteratorsMD

import Base: start, _start, done, next, getindex, setindex!, linearindexing
import Base: @nref, @ncall, @nif, @nexprs, LinearFast, LinearSlow

export eachindex

# Traits for linear indexing
linearindexing(::BitArray) = LinearFast()

Jutho

import Base: start, done, next, getindex, setindex!
import Base: @nref, @ncall, @nif, @nexprs

export eachelement, eachindex, linearindexing, LinearFast

# Traits for linear indexing
abstract LinearIndexing
immutable LinearFast <: LinearIndexing end
immutable LinearSlow <: LinearIndexing end

Jutho

using Base.Cartesian
function permutedimsnew!{T1,T2,N}(P::StridedArray{T1,N},B::StridedArray{T2,N},perm)
    length(perm) == N || error("expected permutation of size $N, but length(perm)=$(length(perm))")
    isperm(perm) || error("input is not a permutation")
    dims = size(P)
    for i = 1:N
        dims[i] == size(B,perm[i]) || throw(DimensionMismatch("destination tensor of incorrect size"))
    end
    stridesB = strides(B)[perm]
    stridesP = strides(P)