From e93997dfac76281ca5f52a27067a567e66c4956f Mon Sep 17 00:00:00 2001
From: nomissbowling <nomissbowling@yahoo.co.jp>
Date: Fri, 21 Oct 2022 16:40:56 +0900
Subject: [PATCH] new branch ComplexLA_develop to solve complex eigenvalues
 eigenvectors and inv etc for general (symmetry and asymmetry) Matrix

  neo/cla.nim
  neo/clasolve.nim
  neo/private/clacommon.nim
  neo/private/clareim.nim
  tests/tcla.nim
  tests/cla/cbase.nim
  tests/cla/cpauli.nim
  tests/cla/chermitian.nim
  tests/cla/cdet.nim
  tests/cla/cinv.nim
  tests/cla/csolve.nim
  tests/cla/ceig.nim

use [c|z]axpy [c|z]dotu [c|z]dotc [c|z]gemv [c|z]gemm [c|z]scal (nimblas)
use [c|z]geev [c|z]gesv overload with macro(czfortran) (nimlapack)
include nimlapack instead of import nimlapack to use private types
template czelColMajor, czelRowMajor for cast[ptr Complex[T]]
unary operator `-`
template compareApprox[T](x, y: complex.Complex[T]): bool

proc all*[T](v: Vector[Complex[T]];
  pred: proc(i: int, x: complex.Complex[T]): bool {.closure.}): bool
proc all*[T](a: Matrix[Complex[T]];
  pred: proc(i, j: int, x: complex.Complex[T]): bool {.closure.}): bool

template `!=~`*[T](x, y: complex.Complex[T]): bool
template `!=~`*[T](v, w: Vector[Complex[T]]): bool
template `!=~`*[T](a, b: Matrix[Complex[T]]): bool

scal for operator `*='*[T](v: var Vector[Complex[T]]; c: complex.Complex[T])
scal for operator `*='*[T](a: var Matrix[Complex[T]]; c: complex.Complex[T])

dotu dotc[Vector[Complex[T]]] Hermitian Inner Product
trace and tr and dotc[Matrix[Complex[T]]] Hilbert-Schmidt Inner Product

TODO: support and check for sparse matrix
TODO: support and check for cuda cudadense cudasparse
---
 neo/cla.nim               | 393 ++++++++++++++++++++++
 neo/clasolve.nim          |  79 +++++
 neo/private/clacommon.nim |  52 +++
 neo/private/clareim.nim   |  65 ++++
 tests/cla/cbase.nim       | 672 ++++++++++++++++++++++++++++++++++++++
 tests/cla/cdet.nim        |  32 ++
 tests/cla/ceig.nim        | 139 ++++++++
 tests/cla/chermitian.nim  |  36 ++
 tests/cla/cinv.nim        |  76 +++++
 tests/cla/cpauli.nim      |  76 +++++
 tests/cla/csolve.nim      |  55 ++++
 tests/tcla.nim            |   9 +
 12 files changed, 1684 insertions(+)
 create mode 100644 neo/cla.nim
 create mode 100644 neo/clasolve.nim
 create mode 100644 neo/private/clacommon.nim
 create mode 100644 neo/private/clareim.nim
 create mode 100644 tests/cla/cbase.nim
 create mode 100644 tests/cla/cdet.nim
 create mode 100644 tests/cla/ceig.nim
 create mode 100644 tests/cla/chermitian.nim
 create mode 100644 tests/cla/cinv.nim
 create mode 100644 tests/cla/cpauli.nim
 create mode 100644 tests/cla/csolve.nim
 create mode 100644 tests/tcla.nim

diff --git a/neo/cla.nim b/neo/cla.nim
new file mode 100644
index 0000000..f795638
--- /dev/null
+++ b/neo/cla.nim
@@ -0,0 +1,393 @@
+# complex linear algebra
+#
+# see also clasolve.nim
+#
+
+import nimblas
+import sequtils, complex
+import sugar
+import ./dense
+
+template czelColMajor[T](ap: CPointer[T]; a, i, j: untyped): untyped=
+  (cast[ptr Complex[T]](ap[2 * (j * a.ld + i)].addr))[] # complex size
+
+template czelRowMajor[T](ap: CPointer[T]; a, i, j: untyped): untyped=
+  (cast[ptr Complex[T]](ap[2 * (i * a.ld + j)].addr))[] # complex size
+
+proc all*[T](v: Vector[Complex[T]];
+  pred: proc(i: int; x: complex.Complex[T]): bool {.closure.}): bool=
+  result = true
+  for i, e in v:
+    if not pred(i, e):
+      result = false
+      break
+
+proc all*[T](a: Matrix[Complex[T]];
+  pred: proc(i, j: int; x: complex.Complex[T]): bool {.closure.}): bool=
+  result = true
+  for ij, e in a:
+    if not pred(ij[0], ij[1], e):
+      result = false
+      break
+
+template tr*[T](a: Matrix[T]): T=
+  a.trace
+
+template tr*[T](a: Matrix[Complex[T]]): complex.Complex[T]=
+  a.trace
+
+template trace*[T](a: Matrix[T]): T=
+  a.vdiag.sum
+
+template trace*[T](a: Matrix[Complex[T]]): complex.Complex[T]=
+  a.vdiag.sum
+
+proc vdiag*[T](a: Matrix[T]): Vector[T] {.noInit.}=
+  let k = min(a.M, a.N)
+  result = zeros(k, T)
+  for i in 0..<k: result[i] = a[i, i]
+
+proc vdiag*[T](a: Matrix[Complex[T]]): Vector[Complex[T]] {.noInit.}=
+  let k = min(a.M, a.N)
+  result = constantVector(k, complex.complex[T](0.0, 0.0))
+  for i in 0..<k: result[i] = a[i, i]
+
+proc diag*[T](v: Vector[T], order=colMajor): Matrix[T] {.noInit.}=
+  result = zeros(v.len, v.len, T, order=order)
+  result.shape.incl(Diagonal)
+  for i in 0..<v.len: result[i, i] = v[i]
+
+proc diag*[T](v: Vector[Complex[T]], order=colMajor): Matrix[Complex[T]] {.noInit.}=
+  result = constantMatrix(v.len, v.len, complex.complex[T](0.0, 0.0), order=order)
+  result.shape.incl(Diagonal)
+  for i in 0..<v.len: result[i, i] = v[i]
+
+proc realize*[T](c: complex.Complex[T]):
+  tuple[real: T, img: T]=
+  result = (c.re, c.im)
+
+proc realize*[T](v: Vector[Complex[T]]):
+  tuple[real: Vector[T], img: Vector[T]]=
+  result = (v.map(x => x.re), v.map(x => x.im))
+
+proc realize*[T](a: Matrix[Complex[T]]):
+  tuple[real: Matrix[T], img: Matrix[T]]=
+  result = (a.map(x => x.re), a.map(x => x.im))
+
+proc toComplex*[T](re: T, im: T): complex.Complex[T]=
+  result = complex.complex[T](re, im)
+
+proc toComplex*[T](re: T): complex.Complex[T]=
+  result = complex.complex[T](re, 0.0)
+
+proc toComplex*[T](re: Vector[T], im: Vector[T]): Vector[Complex[T]]=
+  result = makeVector(re.len,
+    proc(i: int): Complex[T]= complex.complex[T](re[i], im[i]))
+
+proc toComplex*[T](re: Vector[T]): Vector[Complex[T]]=
+  result = re.toComplex(zeros(re.len, T))
+
+proc toComplex*[T](re: Matrix[T], im: Matrix[T]): Matrix[Complex[T]]=
+  result = makeMatrix(re.M, re.N,
+    proc(i, j: int): Complex[T]= complex.complex[T](re[i, j], im[i, j]),
+    order=re.order)
+
+proc toComplex*[T](re: Matrix[T]): Matrix[Complex[T]]=
+  result = re.toComplex(zeros(re.M, re.N, T, order=re.order))
+
+proc to32*(c: complex.Complex[float64]): complex.Complex[float32]=
+  let (re, im) = c.realize
+  result = toComplex[float32](re, im) # not use re.toComplex[float32](im)
+
+proc to64*(c: complex.Complex[float32]): complex.Complex[float64]=
+  let (re, im) = c.realize
+  result = toComplex[float64](re, im) # not use re.toComplex[float64](im)
+
+proc to32*(v: Vector[Complex[float64]]): Vector[Complex[float32]]=
+  let (re, im) = v.realize
+  result = re.to32.toComplex(im.to32)
+
+proc to64*(v: Vector[Complex[float32]]): Vector[Complex[float64]]=
+  let (re, im) = v.realize
+  result = re.to64.toComplex(im.to64)
+
+proc to32*(a: Matrix[Complex[float64]]): Matrix[Complex[float32]]=
+  let (re, im) = a.realize
+  result = re.to32.toComplex(im.to32)
+
+proc to64*(a: Matrix[Complex[float32]]): Matrix[Complex[float64]]=
+  let (re, im) = a.realize
+  result = re.to64.toComplex(im.to64)
+
+proc conjugate*[T](v: Vector[Complex[T]]): Vector[Complex[T]]=
+  result = v.map(x => x.conjugate)
+
+proc conjugate*[T](a: Matrix[Complex[T]]): Matrix[Complex[T]]=
+  result = a.map(x => x.conjugate)
+
+proc H*[T](a: Matrix[Complex[T]]): Matrix[Complex[T]]=
+  result = a.T.conjugate # or a.conjugate.T
+
+template compareApprox[T](x, y: complex.Complex[T]): bool=
+  # TODO: more fast
+  const epsilon = 0.000001
+  const czepsilonabs2 = complex.complex[T](epsilon, epsilon).abs2
+  (x - y).abs2 < czepsilonabs2
+
+proc `=~`*[T](x, y: complex.Complex[T]): bool=
+  result = compareApprox(x, y)
+
+template `!=~`*[T](x, y: complex.Complex[T]): bool=
+  not (x =~ y)
+
+proc `=~`*[T](v, w: Vector[Complex[T]]): bool=
+  # TODO: more fast
+  result = v.all(proc(i: int; x: complex.Complex[T]): bool= x =~ w[i])
+
+template `!=~`*[T](v, w: Vector[Complex[T]]): bool=
+  not (v =~ w)
+
+proc `=~`*[T](a, b: Matrix[Complex[T]]): bool=
+  # TODO: more fast
+  result = a.all(proc(i, j: int; x: complex.Complex[T]): bool= x =~ b[i, j])
+
+template `!=~`*[T](a, b: Matrix[Complex[T]]): bool=
+  not (a =~ b)
+
+proc `+=`*[T](v: var Vector[Complex[T]]; w: Vector[Complex[T]])=
+  assert v.len == w.len
+  var
+    k = v.len
+    alpha = complex.complex[T](1.0, 0.0) # v := w + v
+    vw = w # needless .clone
+    incx = 1
+    incy = 1
+  axpy(k, alpha.addr, vw.fp, incx, v.fp, incy)
+
+proc `+`*[T](v, w: Vector[Complex[T]]): Vector[Complex[T]]=
+  result = v.clone
+  result += w
+
+proc `+=`*[T](a: var  Matrix[Complex[T]]; b: Matrix[Complex[T]])=
+  assert a.M == b.M and a.N == b.N # plus for same shape Matrix
+  if a.isFull and b.isFull and a.order == b.order:
+    var
+      k = a.M * a.N
+      alpha = complex.complex[T](1.0, 0.0) # a := b + a
+      vw = b # needless .clone
+      incx = 1
+      incy = 1
+    axpy(k, alpha.addr, vw.fp, incx, a.fp, incy)
+  else:
+    let ap = cast[CPointer[T]](a.fp) # complex size
+    if a.order == colMajor:
+      for t, x in b:
+        let (i, j) = t
+        czelColMajor(ap, a, i, j) += x
+    else:
+      for t, x in b:
+        let (i, j) = t
+        czelRowMajor(ap, a, i, j) += x
+
+proc `+`*[T](a, b: Matrix[Complex[T]]): Matrix[Complex[T]]=
+  result = a.clone
+  result += b
+
+proc `-=`*[T](v: var Vector[Complex[T]]; w: Vector[Complex[T]])=
+  assert v.len == w.len
+  var
+    k = v.len
+    alpha = complex.complex[T](-1.0, 0.0) # v := -w + v
+    vw = w # needless .clone
+    incx = 1
+    incy = 1
+  axpy(k, alpha.addr, vw.fp, incx, v.fp, incy)
+
+proc `-`*[T](v, w: Vector[Complex[T]]): Vector[Complex[T]]=
+  result = v.clone
+  result -= w
+
+proc `-=`*[T](a: var  Matrix[Complex[T]]; b: Matrix[Complex[T]])=
+  assert a.M == b.M and a.N == b.N # plus for same shape Matrix
+  if a.isFull and b.isFull and a.order == b.order:
+    var
+      k = a.M * a.N
+      alpha = complex.complex[T](-1.0, 0.0) # a := -b + a
+      vw = b # needless .clone
+      incx = 1
+      incy = 1
+    axpy(k, alpha.addr, vw.fp, incx, a.fp, incy)
+  else:
+    let ap = cast[CPointer[T]](a.fp) # complex size
+    if a.order == colMajor:
+      for t, x in b:
+        let (i, j) = t
+        czelColMajor(ap, a, i, j) -= x
+    else:
+      for t, x in b:
+        let (i, j) = t
+        czelRowMajor(ap, a, i, j) -= x
+
+proc `-`*[T](a, b: Matrix[Complex[T]]): Matrix[Complex[T]]=
+  result = a.clone
+  result -= b
+
+proc `*`*[T](v: Vector[T]; c: complex.Complex[T]): Vector[Complex[T]]=
+  result = v.toComplex
+  result *= c
+
+proc `*=`*[T](v: var Vector[Complex[T]]; c: complex.Complex[T])=
+  var
+    k = v.len
+    alpha = c # v := cv
+    incx = 1
+  scal(k, alpha.addr, v.fp, incx)
+
+proc `*`*[T](v: Vector[Complex[T]]; c: complex.Complex[T]):
+  Vector[Complex[T]]=
+  result = v.clone
+  result *= c
+
+proc `*`*[T](v, w: Vector[Complex[T]]): complex.Complex[T]=
+  result = v.dotu(w)
+
+proc dotu*[T](v, w: Vector[Complex[T]]): complex.Complex[T]=
+  assert v.len == w.len
+  var
+    k = v.len
+    vv = v # needless .clone
+    vw = w # needless .clone
+    incx = 1
+    incy = 1
+  result = complex.complex[T](0.0, 0.0)
+  dotu(k, vv.fp, incx, vw.fp, incy, result.addr)
+
+proc dotc*[T](v, w: Vector[Complex[T]]): complex.Complex[T]=
+  # v.dotc(w) == v.conjugate.dotu(w)
+  assert v.len == w.len
+  var
+    k = v.len
+    vv = v # needless .clone
+    vw = w # needless .clone
+    incx = 1
+    incy = 1
+  result = complex.complex[T](0.0, 0.0)
+  dotc(k, vv.fp, incx, vw.fp, incy, result.addr)
+
+proc `*`*[T](a: Matrix[T]; c: complex.Complex[T]): Matrix[Complex[T]]=
+  result = a.toComplex
+  result *= c
+
+proc `*=`*[T](a: var Matrix[Complex[T]]; c: complex.Complex[T])=
+  var
+    k = a.M * a.N
+    alpha = c # a := ca
+    incx = 1
+  scal(k, alpha.addr, a.fp, incx)
+
+proc `*`*[T](a: Matrix[Complex[T]]; c: complex.Complex[T]):
+  Matrix[Complex[T]]=
+  result = a.clone
+  result *= c
+
+proc `*=`*[T](a: var Matrix[Complex[T]]; b: Matrix[Complex[T]])=
+  a = a * b
+
+proc `*`*[T](a, b: Matrix[Complex[T]]): Matrix[Complex[T]]=
+  assert a.N == b.M
+  var
+    m = a.M
+    n = b.N
+    k = b.M # == a.N
+    alpha = complex.complex[T](1.0, 0.0)
+    va = a # needless .clone
+    vb = b # needless .clone
+    beta = complex.complex[T](-1.0, 0.0)
+  result = constantMatrix(m, n, complex.complex[T](0.0, 0.0), order=a.order)
+  if a.order == colMajor and b.order == colMajor:
+    gemm(result.order, noTranspose, noTranspose, m, n, k,
+      alpha.addr, va.fp, va.ld, vb.fp, vb.ld,
+      beta.addr, result.fp, result.ld) # result colMajor
+  elif a.order == rowMajor and b.order == rowMajor:
+    gemm(result.order, noTranspose, noTranspose, m, n, k,
+      alpha.addr, va.fp, va.ld, vb.fp, vb.ld,
+      beta.addr, result.fp, result.ld) # result rowMajor
+  elif a.order == colMajor and b.order == rowMajor:
+    gemm(result.order, noTranspose, transpose, m, n, k,
+      alpha.addr, va.fp, va.ld, vb.fp, vb.ld,
+      beta.addr, result.fp, result.ld) # result colMajor
+  else: # a.order == rowMajor and b.order == colMajor
+    result.order = colMajor
+    result.ld = m
+    gemm(result.order, transpose, noTranspose, m, n, k,
+      alpha.addr, va.fp, va.ld, vb.fp, vb.ld,
+      beta.addr, result.fp, result.ld) # result colMajor != a.order
+
+template dotc*[T](a, b: Matrix[Complex[T]]): complex.Complex[T]=
+  # Hilbert-Schmidt Inner Product
+  result = (a.H * b).tr
+
+proc `*`*[T](a: Matrix[Complex[T]]; v: Vector[Complex[T]]): Vector[Complex[T]]=
+  assert a.N == v.len
+  var
+    m = a.M
+    n = v.len # == a.N
+    alpha = complex.complex[T](1.0, 0.0)
+    va = a # needless .clone
+    vv = v # needless .clone
+    beta = complex.complex[T](-1.0, 0.0)
+    incx = 1
+    incy = 1
+  result = constantVector(n, complex.complex[T](0.0, 0.0))
+  gemv(a.order, noTranspose, m, n,
+    alpha.addr, va.fp, m, vv.fp, incx,
+    beta.addr, result.fp, incy)
+
+template `*`*[T](c: complex.Complex[T],
+  v: Vector[T] or Matrix[T]): auto= v * c
+
+template `*`*[T](c: complex.Complex[T],
+  v: Vector[Complex[T]] or Matrix[Complex[T]]): auto= v * c
+
+template `*=`*[T](v: var Vector[Complex[T]] or var Matrix[Complex[T]]; k: T)=
+  v *= k.toComplex
+
+template `*`*[T](v: Vector[Complex[T]] or Matrix[Complex[T]]; k: T): auto=
+  v * k.toComplex
+
+template `*`*[T](k: T,
+  v: Vector[Complex[T]] or Matrix[Complex[T]]): auto= v * k
+
+template `/=`*[T](v: var Vector[T] or Matrix[T],
+  c: complex.Complex[T])=
+  v *= (complex.complex[T](1.0, 0.0) / c)
+
+template `/`*[T](v: Vector[T] or Matrix[T],
+  c: complex.Complex[T]): auto=
+  v * (complex.complex[T](1.0, 0.0) / c)
+
+template `/=`*[T](v: var Vector[Complex[T]] or Matrix[Complex[T]],
+  c: complex.Complex[T])=
+  v *= (complex.complex[T](1.0, 0.0) / c)
+
+template `/`*[T](v: Vector[Complex[T]] or Matrix[Complex[T]],
+  c: complex.Complex[T]): auto=
+  v * (complex.complex[T](1.0, 0.0) / c)
+
+template `/=`*[T](v: var Vector[Complex[T]] or var Matrix[Complex[T]]; k: T)=
+  v /= k.toComplex
+
+template `/`*[T](v: Vector[Complex[T]] or Matrix[Complex[T]]; k: T): auto=
+  v / k.toComplex
+
+template sameVector*[T](a, b: Matrix[T]): bool=
+  a.asVector == b.asVector
+
+template sameVector*[T](a, b: Matrix[Complex[T]]): bool=
+  a.asVector == b.asVector
+
+template `-`*[T](v: Vector[Complex[T]] or Matrix[Complex[T]]): auto=
+  v * complex.complex[T](-1.0, 0.0)
+
+template `-`*[T](v: Vector[T] or Matrix[T]): auto=
+  v * -1.0
diff --git a/neo/clasolve.nim b/neo/clasolve.nim
new file mode 100644
index 0000000..0a99ada
--- /dev/null
+++ b/neo/clasolve.nim
@@ -0,0 +1,79 @@
+# complex linear algebra solvers
+#
+# see also cla.nim
+#
+
+import sequtils, complex
+import ./dense, ./private/neocommon
+import ./cla, ./private/clacommon
+
+## use czfortran() for lapack_complex_* forked from fortran() in neocommon.nim
+overload(geev, cgeev, zgeev)
+overload(gesv, cgesv, zgesv)
+
+proc usvd*[T: SomeFloat](a: Matrix[Complex[T]]):
+  tuple[U: Matrix[Complex[T]], S: Matrix[Complex[T]], Vh: Matrix[Complex[T]]]=
+  var
+    m = a.M.int32
+    n = a.N.int32
+    vS = a.clone # overwritten by result (Diag or Triangular)
+    vW = constantVector(n, complex.complex[T](0.0, 0.0)) # W: EigenValues
+    jobvl = cstring"V"
+    ldvl = m
+    vVL = constantMatrix(m, m, complex.complex[T](0.0, 0.0), order=a.order)
+    jobvr = cstring"V"
+    ldvr = n
+    vVR = constantMatrix(n, n, complex.complex[T](0.0, 0.0), order=a.order) # VR: EigenVectors
+    lwork = n * 4 # or get before 2pass [c|z]geev
+    vwork = constantVector(lwork, complex.complex[T](0.0, 0.0))
+    rwork = zeros(lwork, T)
+    info = 0'i32
+  czfortran(geev, jobvl, jobvr, n, vS, m, vW,
+    vVL, ldvl, vVR, ldvr, vwork, lwork, rwork, info)
+  # TODO: convert vS to vS.vdiag and vVL to vVR.inv for compatible GE SVD
+  assert vW == vS.vdiag
+  result = (vVL, vS, vVR)
+
+proc usvd*[T: SomeFloat](a: Matrix[T]):
+  tuple[U: Matrix[Complex[T]], S: Matrix[Complex[T]], Vh: Matrix[Complex[T]]]=
+  result = a.toComplex.usvd
+
+proc geneig*[T](a: Matrix[Complex[T]]): (EigenValues[T], EigenVectors[T])=
+  let (_, S, Vh) = a.usvd # (U, S, Vh)
+  var (ers, eis) = (newSeq[T](a.N), newSeq[T](a.N))
+  var (vrs, vis) = (newSeq[Vector[T]](a.N), newSeq[Vector[T]](a.N))
+  for i in 0..<a.N:
+    (ers[i], eis[i]) = S[i, i].realize
+    (vrs[i], vis[i]) = Vh.column(i).realize # .T
+  result = (
+    EigenValues[T](real: ers, img: eis),
+    EigenVectors[T](real: vrs, img: vis))
+
+proc geneig*[A: SomeFloat](a: Matrix[A]): (EigenValues[A], EigenVectors[A])=
+  result = a.toComplex.geneig
+
+proc solve*[T: SomeFloat](a, b: Matrix[Complex[T]]): Matrix[Complex[T]]=
+  var
+    m = a.M.int32
+    n = a.N.int32 # == b.M
+    nrhs = b.N.int32
+    va = a.clone
+    ipiv = newSeq[int32](n)
+    info = 0'i32
+  result = b.clone
+  czfortran(gesv, n, nrhs, va, m, ipiv, result, n, info)
+
+proc solve*[T: SomeFloat](a: Matrix[Complex[T]]; v: Vector[Complex[T]]):
+  Vector[Complex[T]]=
+  result = a.solve(matrix(@[v.toSeq], order=a.order).t).column(0)
+
+proc inv*[T](a: Matrix[Complex[T]]): Matrix[Complex[T]]=
+  result = a.solve(eye(a.N, T, order=a.order).toComplex)
+
+proc det*[T](a: Matrix[Complex[T]]): complex.Complex[T]=
+  # assert a.M == a.N # TODO: must be assertion error ?
+  if a.M != a.N: result = complex.complex[T](0.0, 0.0)
+  else:
+    let (_, S, _) = a.usvd # (U, S, Vh)
+    result = complex.complex[T](1.0, 0.0)
+    for i in 0..<a.N: result *= S[i, i]
diff --git a/neo/private/clacommon.nim b/neo/private/clacommon.nim
new file mode 100644
index 0000000..9f38071
--- /dev/null
+++ b/neo/private/clacommon.nim
@@ -0,0 +1,52 @@
+# complex linear algebra solvers
+
+import macros
+import complex
+include nimlapack # (not use import) private type lapack_complex_[float|double]
+
+## cast for [c|z]geev [c|z]gesv in nimlapack.nim
+## to fix .fp ambiguous between ptr lapack_complex_[float|double]
+##  Vector or Matrix[Complex[float32 | float64]]
+
+template force_cast_ptr*[P: complex.Complex[float32]](p: ptr P):
+  ptr lapack_complex_float=
+  cast[ptr lapack_complex_float](p)
+
+template force_cast_ptr*[P: complex.Complex[float64]](p: ptr P):
+  ptr lapack_complex_double=
+  cast[ptr lapack_complex_double](p)
+
+template force_cast_ptr*[P: float32](p: ptr P):
+  ptr cfloat=
+  cast[ptr cfloat](p)
+
+template force_cast_ptr*[P: float64](p: ptr P):
+  ptr cdouble=
+  cast[ptr cdouble](p)
+
+template force_cast_ptr*[T, P](p: ptr P): ptr T=
+  cast[ptr T](p)
+
+proc czgetAddress(n: NimNode): NimNode =
+  let t = n.getTypeImpl
+  if t.kind == nnkBracketExpr and $(t[0]) == "seq":
+    result = quote do:
+      addr `n`[0]
+  elif t.kind == nnkRefTy and $(t[0]) in ["Vector", "Vector:ObjectType"]:
+    result = quote do:
+      force_cast_ptr(`n`.fp)
+  elif t.kind == nnkRefTy and $(t[0]) in ["Matrix", "Matrix:ObjectType"]:
+    result = quote do:
+      force_cast_ptr(`n`.fp)
+  elif $t == "cstring":
+    result = quote do:
+      `n`
+  else:
+    result = quote do:
+      addr `n`
+
+macro czfortran*(f: untyped, callArgs: varargs[typed]): auto =
+  var transformedCallArgs = newSeqOfCap[NimNode](callArgs.len)
+  for x in callArgs:
+    transformedCallArgs.add(czgetAddress(x))
+  result = newCall(f, transformedCallArgs)
diff --git a/neo/private/clareim.nim b/neo/private/clareim.nim
new file mode 100644
index 0000000..eadcf13
--- /dev/null
+++ b/neo/private/clareim.nim
@@ -0,0 +1,65 @@
+# complex linear algebra privates (re im decomposition)
+
+import complex, ../dense, ../cla
+
+proc ae[T](x, y: complex.Complex[T]): bool=
+  result = x =~ y
+
+proc ae[T](v, w: Vector[Complex[T]]): bool=
+  result = v =~ w
+
+proc ae_ri[T](v, w: Vector[Complex[T]]): bool=
+  let (vr, vi) = v.realize
+  let (wr, wi) = w.realize
+  result = vr =~ wr and vi =~ wi
+
+proc ae[T](a, b: Matrix[Complex[T]]): bool=
+  result = a =~ b
+
+proc ae_ri[T](a, b: Matrix[Complex[T]]): bool=
+  let (ar, ai) = a.realize
+  let (br, bi) = b.realize
+  result = ar =~ br and ai =~ bi
+
+proc plus[T](v, w: Vector[Complex[T]]): Vector[Complex[T]]=
+  result = v + w
+
+proc plus_ri[T](v, w: Vector[Complex[T]]): Vector[Complex[T]]=
+  let (vr, vi) = v.realize
+  let (wr, wi) = w.realize
+  result = (vr + wr).toComplex(vi + wi)
+
+proc plus[T](a, b: Matrix[Complex[T]]): Matrix[Complex[T]]=
+  result = a + b
+
+proc plus_ri[T](a, b: Matrix[Complex[T]]): Matrix[Complex[T]]=
+  let (ar, ai) = a.realize
+  let (br, bi) = b.realize
+  result = (ar + br).toComplex(ai + bi)
+
+proc dotu_ri[T](v, w: Vector[Complex[T]]): complex.Complex[T]=
+  let (vr, vi) = v.realize
+  let (wr, wi) = w.realize
+  result = (vr * wr - vi * wi).toComplex(vr * wi + vi * wr)
+
+proc dotc_ri[T](v, w: Vector[Complex[T]]): complex.Complex[T]=
+  let (vr, vi) = v.realize
+  let (wr, wi) = w.realize
+  result = (vr * wr + vi * wi).toComplex(vr * wi - vi * wr)
+
+proc dotu[T](a, b: Matrix[Complex[T]]): Matrix[Complex[T]]=
+  result = a * b
+
+proc dotu_ri[T](a, b: Matrix[Complex[T]]): Matrix[Complex[T]]=
+  let (ar, ai) = a.realize
+  let (br, bi) = b.realize
+  result = (ar * br - ai * bi).toComplex(ar * bi + ai * br)
+
+proc dotu[T](a: Matrix[Complex[T]]; v: Vector[Complex[T]]): Vector[Complex[T]]=
+  result = a * v
+
+proc dotu_ri[T](a: Matrix[Complex[T]]; v: Vector[Complex[T]]):
+  Vector[Complex[T]]=
+  let (ar, ai) = a.realize
+  let (vr, vi) = v.realize
+  result = (ar * vr - ai * vi).toComplex(ar * vi + ai * vr)
diff --git a/tests/cla/cbase.nim b/tests/cla/cbase.nim
new file mode 100644
index 0000000..7c94e44
--- /dev/null
+++ b/tests/cla/cbase.nim
@@ -0,0 +1,672 @@
+# complex linear algebra solvers
+
+import unittest, complex, neo/cla, neo/dense
+
+proc run() =
+  suite "basic complex computations":
+    # TODO: assert
+    # TODO: support and check for sparse matrix
+    # TODO: support and check for cuda cudadense cudasparse
+
+    const
+      ci = im(1.0)
+      cmi = im(-1.0)
+      coi = complex.complex(1.0, 1.0)
+      cmomi = complex.complex(-1.0, -1.0)
+      ci32 = ci.to32
+      cmi32 = cmi.to32
+      coi32 = coi.to32
+      cmomi32 = cmomi.to32
+
+    test "Complex[float32] Complex[float64] of a complex vector matrix":
+      let
+        vf32re = vector(@[0.1'f32, 0.2'f32])
+        vf32 = vf32re.toComplex
+        mf32re = matrix(@[
+          @[0.1'f32, 0.4'f32],
+          @[0.2'f32, 0.3'f32]
+        ])
+        mf32 = mf32re.toComplex
+        vf64re = vector(@[0.1'f64, 0.2'f64])
+        vf64 = vf64re.toComplex
+        mf64re = matrix(@[
+          @[0.1'f64, 0.4'f64],
+          @[0.2'f64, 0.3'f64]
+        ])
+        mf64 = mf64re.toComplex
+        vf64delta5 = vector(@[0.10001'f64, 0.19999'f64]).toComplex
+        mf64delta5 = matrix(@[
+          @[0.10001'f64, 0.39999'f64],
+          @[0.19999'f64, 0.30001'f64]
+        ]).toComplex
+        vf64delta6 = vector(@[0.100001'f64, 0.199999'f64]).toComplex
+        mf64delta6 = matrix(@[
+          @[0.100001'f64, 0.399999'f64],
+          @[0.199999'f64, 0.300001'f64]
+        ]).toComplex
+
+      test "to32 of a complex vector matrix (float64)":
+        # check(vf32 =~ vf64) # type mismatch on compile
+        # check(mf32 =~ mf64) # type mismatch on compile
+        check(vf32 =~ vf64.to32)
+        check(mf32 =~ mf64.to32)
+
+        check((vf64 * 2.0) == (vf64re * coi * coi.conjugate))
+        check((mf64 * 2.0) == (mf64re * coi * coi.conjugate))
+
+      test "to64 of a complex vector matrix (float32)":
+        check(vf32.to64 =~ vf64)
+        check(mf32.to64 =~ mf64)
+
+        check((vf32 * 2.0'f32) == (vf32re * coi32 * coi32.conjugate))
+        check((mf32 * 2.0'f32) == (mf32re * coi32 * coi32.conjugate))
+
+      test "epsilon precision of a complex vector matrix (float64)":
+        check(not (vf64delta5 == vf64))
+        check(not (vf64delta5 =~ vf64)) # epsilon precision
+        check(not (mf64delta5 == mf64))
+        check(not (mf64delta5 =~ mf64)) # epsilon precision
+
+        check(not (vf64delta6 == vf64))
+        check(vf64delta6 =~ vf64) # epsilon precision
+        check(not (mf64delta6 == mf64))
+        check(mf64delta6 =~ mf64) # epsilon precision
+
+      test "epsilon precision of a complex vector matrix (float32)":
+        let
+          vf32delta5 = vf64delta5.to32
+          mf32delta5 = mf64delta5.to32
+          vf32delta6 = vf64delta6.to32
+          mf32delta6 = mf64delta6.to32
+
+        check(not (vf32delta5 == vf32))
+        check(not (vf32delta5 =~ vf32)) # epsilon precision
+        check(not (mf32delta5 == mf32))
+        check(not (mf32delta5 =~ mf32)) # epsilon precision
+
+        check(not (vf32delta6 == vf32))
+        check(vf32delta6 =~ vf32) # epsilon precision
+        check(not (mf32delta6 == mf32))
+        check(mf32delta6 =~ mf32) # epsilon precision
+
+    test "unary operator or real times of a complex vector matrix":
+      let
+        v3are = vector(@[-1.0, 0.1, -0.1])
+        v3a = v3are.toComplex
+        v3mare = vector(@[1.0, -0.1, 0.1])
+        v3ma = v3mare.toComplex
+        m22are = matrix(@[@[0.0, -1.0], @[1.0, 0.0]])
+        m22a = m22are.toComplex
+        m22mare = matrix(@[@[0.0, 1.0], @[-1.0, 0.0]])
+        m22ma = m22mare.toComplex
+
+      test "unary operator or real times of a complex vector matrix (float64)":
+        check(0.0 - coi == cmomi) # OK only for scalar
+        check((-coi) == cmomi)
+        check(-1.0 * coi == cmomi)
+
+        ## type mismatch (at position: 2)
+        # check(0.0 - v3are == v3mare)
+        # check(0.0 - v3a == v3ma)
+        # check(0.0 - m22are == m22mare)
+        # check(0.0 - m22a == m22ma)
+
+        check((-v3are) == v3mare)
+        check((-v3a) == v3ma)
+        check((-m22are) == m22mare)
+        check((-m22a) == m22ma)
+
+        check(-1.0 * v3are == v3mare)
+        check(-1.0 * v3a == v3ma)
+        check(-1.0 * m22are == m22mare)
+        check(-1.0 * m22a == m22ma)
+
+        check(v3are / -1.0 == v3mare)
+        check(v3a / -1.0 == v3ma)
+        check(m22are / -1.0 == m22mare)
+        check(m22a / -1.0 == m22ma)
+
+        var
+          v3vre = v3are.clone
+          v3v = v3a.clone
+          m22vre = m22are.clone
+          m22v = m22a.clone
+          vd3vre = v3are.clone
+          vd3v = v3a.clone
+          md22vre = m22are.clone
+          md22v = m22a.clone
+        v3vre *= -1.0
+        v3v *= -1.0
+        m22vre *= -1.0
+        m22v *= -1.0
+        vd3vre *= -1.0
+        vd3v *= -1.0
+        md22vre *= -1.0
+        md22v *= -1.0
+        check(v3vre == v3mare)
+        check(v3v == v3ma)
+        check(m22vre == m22mare)
+        check(m22v == m22ma)
+        check(vd3vre == v3mare)
+        check(vd3v == v3ma)
+        check(md22vre == m22mare)
+        check(md22v == m22ma)
+
+      test "unary operator or real times of a complex vector matrix (float32)":
+        let
+          v3are32 = v3are.to32
+          v3a32 = v3a.to32
+          v3mare32 = v3mare.to32
+          v3ma32 = v3ma.to32
+          m22are32 = m22are.to32
+          m22a32 = m22a.to32
+          m22mare32 = m22mare.to32
+          m22ma32 = m22ma.to32
+
+        check(0.0'f32 - coi32 == cmomi32) # OK only for scalar
+        check((-coi32) == cmomi32)
+        check(-1.0'f32 * coi32 == cmomi32)
+
+        ## type mismatch (at position: 2)
+        # check(0.0'f32 - v3are32 == v3mare32)
+        # check(0.0'f32 - v3a32 == v3ma32)
+        # check(0.0'f32 - m22are32 == m22mare32)
+        # check(0.0'f32 - m22a32 == m22ma32)
+
+        check((-v3are32) == v3mare32)
+        check((-v3a32) == v3ma32)
+        check((-m22are32) == m22mare32)
+        check((-m22a32) == m22ma32)
+
+        check(-1.0'f32 * v3are32 == v3mare32)
+        check(-1.0'f32 * v3a32 == v3ma32)
+        check(-1.0'f32 * m22are32 == m22mare32)
+        check(-1.0'f32 * m22a32 == m22ma32)
+
+        check(v3are32 / -1.0'f32 == v3mare32)
+        check(v3a32 / -1.0'f32 == v3ma32)
+        check(m22are32 / -1.0'f32 == m22mare32)
+        check(m22a32 / -1.0'f32 == m22ma32)
+
+        var
+          v3vre32 = v3are32.clone
+          v3v32 = v3a32.clone
+          m22vre32 = m22are32.clone
+          m22v32 = m22a32.clone
+          vd3vre32 = v3are32.clone
+          vd3v32 = v3a32.clone
+          md22vre32 = m22are32.clone
+          md22v32 = m22a32.clone
+        v3vre32 *= -1.0'f32
+        v3v32 *= -1.0'f32
+        m22vre32 *= -1.0'f32
+        m22v32 *= -1.0'f32
+        vd3vre32 /= -1.0'f32
+        vd3v32 /= -1.0'f32
+        md22vre32 /= -1.0'f32
+        md22v32 /= -1.0'f32
+        check(v3vre32 == v3mare32)
+        check(v3v32 == v3ma32)
+        check(m22vre32 == m22mare32)
+        check(m22v32 == m22ma32)
+        check(vd3vre32 == v3mare32)
+        check(vd3v32 == v3ma32)
+        check(md22vre32 == m22mare32)
+        check(md22v32 == m22ma32)
+
+    test "order colMajor rowMajor":
+      let
+        v22mem = vector(@[1.0, 2.0, 3.0, 4.0]).toComplex
+        m22default = v22mem.clone.asMatrix(2, 2) # 13|24
+        m22c = v22mem.clone.asMatrix(2, 2, order=colMajor) # 13|24
+        m22r = v22mem.clone.asMatrix(2, 2, order=rowMajor) # 12|34
+
+        v6mem23r = vector(@[1.0, 3.0, 5.0, 2.0, 4.0, 6.0]).toComplex
+        v6mem32r = vector(@[1.0, 4.0, 2.0, 5.0, 3.0, 6.0]).toComplex
+        v6mem = vector(@[1.0, 2.0, 3.0, 4.0, 5.0, 6.0]).toComplex
+        n23c = matrix(@[ # 123456
+          @[1.0, 3.0, 5.0],
+          @[2.0, 4.0, 6.0]
+        ], order=colMajor).toComplex
+        n23r = matrix(@[ # 135246
+          @[1.0, 3.0, 5.0],
+          @[2.0, 4.0, 6.0]
+        ], order=rowMajor).toComplex
+        n32c = matrix(@[ # 123456
+          @[1.0, 4.0],
+          @[2.0, 5.0],
+          @[3.0, 6.0]
+        ], order=colMajor).toComplex
+        n32r = matrix(@[ # 142536
+          @[1.0, 4.0],
+          @[2.0, 5.0],
+          @[3.0, 6.0]
+        ], order=rowMajor).toComplex
+
+      test "order colMajor rowMajor (float64)":
+        check(m22c == m22default)
+        check(not (m22c == m22r))
+        check(not (m22c =~ m22r))
+        check(m22c == m22r.t)
+        check(m22c =~ m22r.t)
+        check(m22c.sameVector(m22r)) # 1234
+        check(m22c.asVector == v22mem) # 1234
+        check(m22r.asVector == v22mem) # 1234
+
+        check(n23c.asVector == v6mem) # 123456
+        check(n23r.asVector == v6mem23r) # 135246
+        check(n32c.asVector == v6mem) # 123456
+        check(n32r.asVector == v6mem32r) # 142536
+
+      test "order colMajor rowMajor (float32)":
+        let
+          v22mem32 = v22mem.to32
+          m22default32 = m22default.to32
+          m22c32 = m22c.to32
+          m22r32 = m22r.to32
+          v6mem23r32 = v6mem23r.to32
+          v6mem32r32 = v6mem32r.to32
+          v6mem32 = v6mem.to32
+          n23c32 = n23c.to32
+          n23r32 = n23r.to32
+          n32c32 = n32c.to32
+          n32r32 = n32r.to32
+
+        check(m22c32 == m22default32)
+        check(not (m22c32 == m22r32))
+        check(not (m22c32 =~ m22r32))
+        check(m22c32 == m22r32.t)
+        check(m22c32 =~ m22r32.t)
+        check(m22c32.sameVector(m22r32)) # 1234
+        check(m22c32.asVector == v22mem32) # 1234
+        check(m22r32.asVector == v22mem32) # 1234
+
+        check(n23c32.asVector == v6mem32) # 123456
+        check(n23r32.asVector == v6mem23r32) # 135246
+        check(n32c32.asVector == v6mem32) # 123456
+        check(n32r32.asVector == v6mem32r32) # 142536
+
+    test "order and diag of a complex matrix":
+      let
+        m23c = matrix(@[
+          @[1.0, 3.0, 5.0],
+          @[2.0, 4.0, 6.0]
+        ], order=colMajor).toComplex
+        m23r = matrix(@[
+          @[1.0, 2.0, 3.0],
+          @[4.0, 5.0, 6.0]
+        ], order=rowMajor).toComplex
+
+        d22c = vector(@[1.0, 2.0]).diag(order=colMajor).toComplex
+        d22r = vector(@[1.0, 2.0]).diag(order=rowMajor).toComplex
+        v23cc = vector(@[1.0, 4.0, 3.0, 8.0, 5.0, 12.0]).toComplex
+        v23cr = vector(@[1.0, 8.0, 2.0, 10.0, 3.0, 12.0]).toComplex
+        v23rc = vector(@[1.0, 4.0, 3.0, 8.0, 5.0, 12.0]).toComplex
+        v23rr = vector(@[1.0, 2.0, 3.0, 8.0, 10.0, 12.0]).toComplex
+
+        d33c = vector(@[1.0, 2.0, 3.0]).diag(order=colMajor).toComplex
+        d33r = vector(@[1.0, 2.0, 3.0]).diag(order=rowMajor).toComplex
+        u23cc = vector(@[1.0, 2.0, 6.0, 8.0, 15.0, 18.0]).toComplex
+        u23rc = vector(@[1.0, 4.0, 4.0, 10.0, 9.0, 18.0]).toComplex
+        u23cr = vector(@[1.0, 2.0, 6.0, 8.0, 15.0, 18.0]).toComplex
+        u23rr = vector(@[1.0, 4.0, 9.0, 4.0, 10.0, 18.0]).toComplex
+
+      test "order and diag of a complex matrix (float64)":
+        check(m23c.sameVector(m23r)) # 123456
+
+        var
+          ms23c = m23c.clone
+          md23c = m23c.clone
+          ms23r = m23r.clone
+          md23r = m23r.clone
+        ms23c *= ci
+        md23c /= cmi
+        ms23r *= ci
+        md23r /= cmi
+        check(ms23c.sameVector(md23c))
+        check(ms23r.sameVector(md23r))
+        check((m23c * ci).sameVector(m23c / cmi))
+        check((m23r * ci).sameVector(m23r / cmi))
+
+        check(d22c == d22r)
+        check((d22c * m23c).asVector == v23cc) # 1 4 3 8 5 12 colMajor
+        check((d22c * m23r).asVector == v23cr) # 1 8 2 10 3 12 colMajor
+        check((d22r * m23c).asVector == v23rc) # 1 4 3 8 5 12 colMajor !
+        check((d22r * m23r).asVector == v23rr) # 1 2 3 8 10 12 rowMajor
+        check((d22c * m23c).sameVector(d22r * m23c)) # 1 4 3 8 5 12
+        check(not (d22c * m23r).sameVector(d22r * m23r))
+
+        check(d33c == d33r)
+        check((m23c * d33c).asVector == u23cc) # 1 2 6 8 15 18 colMajor
+        check((m23r * d33c).asVector == u23rc) # 1 4 4 10 9 18 colMajor !
+        check((m23c * d33r).asVector == u23cr) # 1 2 6 8 15 18 colMajor
+        check((m23r * d33r).asVector == u23rr) # 1 4 9 4 10 18 rowMajor
+        check((m23c * d33c).sameVector(m23c * d33r)) # 1 2 6 8 15 18
+        check(not (m23r * d33c).sameVector(m23r * d33r))
+
+      test "order and diag of a complex matrix (float32)":
+        let
+          m23c32 = m23c.to32
+          m23r32 = m23r.to32
+          d22c32 = d22c.to32
+          d22r32 = d22r.to32
+          v23cc32 = v23cc.to32
+          v23cr32 = v23cr.to32
+          v23rc32 = v23rc.to32
+          v23rr32 = v23rr.to32
+          d33c32 = d33c.to32
+          d33r32 = d33r.to32
+          u23cc32 = u23cc.to32
+          u23rc32 = u23rc.to32
+          u23cr32 = u23cr.to32
+          u23rr32 = u23rr.to32
+
+        check(m23c32.sameVector(m23r32)) # 123456
+
+        var
+          ms23c32 = m23c32.clone
+          md23c32 = m23c32.clone
+          ms23r32 = m23r32.clone
+          md23r32 = m23r32.clone
+        ms23c32 *= ci32
+        md23c32 /= cmi32
+        ms23r32 *= ci32
+        md23r32 /= cmi32
+        check(ms23c32.sameVector(md23c32))
+        check(ms23r32.sameVector(md23r32))
+        check((m23c32 * ci32).sameVector(m23c32 / cmi32))
+        check((m23r32 * ci32).sameVector(m23r32 / cmi32))
+
+        check(d22c32 == d22r32)
+        check((d22c32 * m23c32).asVector == v23cc32) # 1 4 3 8 5 12 colMajor
+        check((d22c32 * m23r32).asVector == v23cr32) # 1 8 2 10 3 12 colMajor
+        check((d22r32 * m23c32).asVector == v23rc32) # 1 4 3 8 5 12 colMajor !
+        check((d22r32 * m23r32).asVector == v23rr32) # 1 2 3 8 10 12 rowMajor
+        check((d22c32 * m23c32).sameVector(d22r32 * m23c32)) # 1 4 3 8 5 12
+        check(not (d22c32 * m23r32).sameVector(d22r32 * m23r32))
+
+        check(d33c32 == d33r32)
+        check((m23c32 * d33c32).asVector == u23cc32) # 1 2 6 8 15 18 colMajor
+        check((m23r32 * d33c32).asVector == u23rc32) # 1 4 4 10 9 18 colMajor !
+        check((m23c32 * d33r32).asVector == u23cr32) # 1 2 6 8 15 18 colMajor
+        check((m23r32 * d33r32).asVector == u23rr32) # 1 4 9 4 10 18 rowMajor
+        check((m23c32 * d33c32).sameVector(m23c32 * d33r32)) # 1 2 6 8 15 18
+        check(not (m23r32 * d33c32).sameVector(m23r32 * d33r32))
+
+    test "vdiag and diag of a complex matrix":
+      let
+        a = matrix(@[
+          @[1.0, 4.0],
+          @[2.0, 3.0]
+        ]).toComplex
+        da = vector(@[1.0, 3.0]).toComplex
+        b = matrix(@[
+          @[0.0, -1.0],
+          @[1.0, 0.0]
+        ]).toComplex
+        db = vector(@[0.0, 0.0]).toComplex
+        m = matrix(@[
+          @[1.0, 0.0, 1.0],
+          @[0.0, 1.0, 2.0]
+        ]).toComplex
+        dm = vector(@[1.0, 1.0]).toComplex
+        n = matrix(@[
+          @[1.0, 0.0],
+          @[0.0, 1.0],
+          @[1.0, 2.0]
+        ]).toComplex
+        dn = vector(@[1.0, 1.0]).toComplex
+
+      test "vdiag and diag of a complex matrix (float64)":
+        check(vdiag(a) =~ da)
+        check(vdiag(b) =~ db)
+        check(vdiag(m) =~ dm) # TODO: must be assertion error ?
+        check(vdiag(n) =~ dn) # TODO: must be assertion error ?
+        check(diag(dm) == diag(dn))
+
+      test "vdiag and diag of a complex matrix (float32)":
+        check(a.to32.vdiag =~ da.to32)
+        check(b.to32.vdiag =~ db.to32)
+        check(m.to32.vdiag =~ dm.to32) # TODO: must be assertion error ?
+        check(n.to32.vdiag =~ dn.to32) # TODO: must be assertion error ?
+        check(dm.to32.diag == dn.to32.diag)
+
+    test "dot scalar and plus minus of a complex vector":
+      let
+        cri = complex.complex(2.0, -1.0)
+        v3cri = vector(@[cri, cri, cri])
+        v3re = vector(@[1.0, -2.0, 3.0])
+        v3im = vector(@[-1.0, 2.0, -3.0])
+        v3 = v3re.toComplex(v3im)
+        w3re = vector(@[-1.0, -2.0, -3.0])
+        w3im = vector(@[1.0, 2.0, -3.0])
+        w3 = w3re.toComplex(w3im)
+        vdc3m = complex.complex(6.0, 2.0)
+        vdc3mreimz = vdc3m.re.toComplex
+        vdu3m = complex.complex(2.0, -6.0)
+        vds3mre = vector(@[1.0, -2.0, 3.0])
+        vds3mim = vector(@[1.0, -2.0, 3.0])
+        vds3m = vds3mre.toComplex(vds3mim)
+        vs3mre = vector(@[1.0, -2.0, 3.0])
+        vs3mim = vector(@[-3.0, 6.0, -9.0])
+        vs3m = vs3mre.toComplex(vs3mim)
+        vpw3mre = vector(@[0.0, -4.0, 0.0])
+        vpw3mim = vector(@[0.0, 4.0, -6.0])
+        vpw3m = vpw3mre.toComplex(vpw3mim)
+        vmw3mre = vector(@[2.0, 0.0, 6.0])
+        vmw3mim = vector(@[-2.0, 0.0, 0.0])
+        vmw3m = vmw3mre.toComplex(vmw3mim)
+
+      test "dot scalar and plus minus of a complex vector (float64)":
+        var
+          vds3v = v3.clone
+          vs3v = v3.clone
+          vpw3v = v3.clone
+          vmw3v = v3.clone
+        vds3v /= cmi
+        vs3v *= cri
+        vpw3v += w3
+        vmw3v -= w3
+        check(vds3v == vds3m)
+        check(vs3v == vs3m)
+        check(vpw3v == vpw3m)
+        check(vmw3v == vmw3m)
+        check(v3.dotc(v3cri) == vdc3m)
+        check(v3.dotc(v3cri) == v3.conjugate.dotu(v3cri))
+        check(v3.dotc(v3cri) == v3cri.dotc(v3).conjugate)
+        check((v3.dotc(v3cri) + v3cri.dotc(v3)) / 2.0 == vdc3mreimz)
+        check(v3.dotu(v3cri) == vdu3m)
+        check(v3 * v3cri == vdu3m)
+        check(v3 * cri == vs3m)
+        check(v3 / cmi == vds3m)
+        check(v3 + w3 == vpw3m)
+        check(v3 - w3 == vmw3m)
+
+      test "dot scalar and plus minus of a complex vector (float32)":
+        let
+          cri32 = cri.to32
+          v3cri32 = v3cri.to32
+          v332 = v3.to32
+          w332 = w3.to32
+          vdc3m32 = vdc3m.to32
+          vdc3mreimz32 = vdc3mreimz.to32
+          vdu3m32 = vdu3m.to32
+          vds3m32 = vds3m.to32
+          vs3m32 = vs3m.to32
+          vpw3m32 = vpw3m.to32
+          vmw3m32 = vmw3m.to32
+        var
+          vds3v32 = v332.clone
+          vs3v32 = v332.clone
+          vpw3v32 = v332.clone
+          vmw3v32 = v332.clone
+        vds3v32 /= cmi32
+        vs3v32 *= cri32
+        vpw3v32 += w332
+        vmw3v32 -= w332
+        check(vds3v32 == vds3m32)
+        check(vs3v32 == vs3m32)
+        check(vpw3v32 == vpw3m32)
+        check(vmw3v32 == vmw3m32)
+        check(v332.dotc(v3cri32) == vdc3m32)
+        check(v332.dotc(v3cri32) == v332.conjugate.dotu(v3cri32))
+        check(v332.dotc(v3cri32) == v3cri32.dotc(v332).conjugate)
+        check((v332.dotc(v3cri32) + v3cri32.dotc(v332)) / 2.0 == vdc3mreimz32)
+        check(v332.dotu(v3cri32) == vdu3m32)
+        check(v332 * v3cri32 == vdu3m32)
+        check(v332 * cri32 == vs3m32)
+        check(v332 / cmi32 == vds3m32)
+        check(v332 + w332 == vpw3m32)
+        check(v332 - w332 == vmw3m32)
+
+    test "order and plus minus of a complex matrix":
+      let
+        p23cim = matrix(@[
+          @[-1.0, -3.0, -5.0],
+          @[-2.0, -4.0, -6.0]
+        ], order=colMajor)
+        p23c = matrix(@[
+          @[1.0, 3.0, 5.0],
+          @[2.0, 4.0, 6.0]
+        ], order=colMajor).toComplex(p23cim)
+        p23rim = matrix(@[
+          @[-1.0, -2.0, -3.0],
+          @[-4.0, -5.0, -6.0]
+        ], order=rowMajor)
+        p23r = matrix(@[
+          @[1.0, 2.0, 3.0],
+          @[4.0, 5.0, 6.0]
+        ], order=rowMajor).toComplex(p23rim)
+        v23mre = vector(@[1.0, 2.0, 3.0, 4.0, 5.0, 6.0])
+        v23mim = vector(@[-1.0, -2.0, -3.0, -4.0, -5.0, -6.0])
+        v23m = v23mre.toComplex(v23mim) # as memory
+        v23cre = vector(@[2.0, 4.0, 6.0, 8.0, 10.0, 12.0])
+        v23cim = vector(@[-2.0, -4.0, -6.0, -8.0, -10.0, -12.0])
+        v23c = v23cre.toComplex(v23cim) # p23r as colMajor
+        v23rre = vector(@[2.0, 6.0, 5.0, 9.0, 8.0, 12.0])
+        v23rim = vector(@[-2.0, -6.0, -5.0, -9.0, -8.0, -12.0])
+        v23r = v23rre.toComplex(v23rim) # p23r as rowMajor
+        v23vre = vector(@[2.0, 5.0, 8.0, 6.0, 9.0, 12.0])
+        v23vim = vector(@[-2.0, -5.0, -8.0, -6.0, -9.0, -12.0])
+        v23v = v23vre.toComplex(v23vim) # p23c as rowMajor
+
+        p23zim = matrix(@[
+          @[-1.0, 2.0, -3.0],
+          @[4.0, -5.0, -6.0]
+        ], order=rowMajor)
+        p23z = matrix(@[
+          @[1.0, -2.0, 3.0],
+          @[-4.0, 5.0, -6.0]
+        ], order=colMajor).toComplex(p23zim)
+        v23cpzre = vector(@[2.0, -2.0, 1.0, 9.0, 8.0, 0.0])
+        v23cpzim = vector(@[-2.0, 2.0, -1.0, -9.0, -8.0, -12.0])
+        v23cpz = v23cpzre.toComplex(v23cpzim)
+        v23rpzre = vector(@[2.0, 0.0, 6.0, 0.0, 10.0, 0.0])
+        v23rpzim = vector(@[-2.0, 0.0, -6.0, 0.0, -10.0, -12.0])
+        v23rpz = v23rpzre.toComplex(v23rpzim)
+        v23cmzre = vector(@[0.0, 6.0, 5.0, -1.0, 2.0, 12.0])
+        v23cmzim = vector(@[0.0, -6.0, -5.0, 1.0, -2.0, 0.0])
+        v23cmz = v23cmzre.toComplex(v23cmzim)
+        v23rmzre = vector(@[0.0, 4.0, 0.0, 8.0, 0.0, 12.0])
+        v23rmzim = vector(@[0.0, -4.0, 0.0, -8.0, 0.0, 0.0])
+        v23rmz = v23rmzre.toComplex(v23rmzim)
+
+      test "order and plus minus of a complex matrix (float64)":
+        check(p23c.asVector == v23m)
+        check(p23r.asVector == v23m)
+        check(not ((p23c + p23r).asVector == v23c))
+        check((p23c + p23r).asVector == v23r)
+        check((p23r + p23c).asVector == v23v)
+        check((p23r + p23r).asVector == v23c)
+        check((p23c + p23c).asVector == v23c)
+
+        var
+          o23cr = p23c.clone
+          o23rc = p23r.clone
+          o23rr = p23r.clone
+          o23cc = p23c.clone
+        o23cr += p23r
+        o23rc += p23c
+        o23rr += p23r
+        o23cc += p23c
+        check(o23cr.asVector == v23r)
+        check(o23rc.asVector == v23v)
+        check(o23rr.asVector == v23c)
+        check(o23cc.asVector == v23c)
+
+        check((p23c + p23z).asVector == v23cpz)
+        check((p23r + p23z).asVector == v23rpz)
+        check((p23c - p23z).asVector == v23cmz)
+        check((p23r - p23z).asVector == v23rmz)
+
+        var
+          o23cpz = p23c.clone
+          o23rpz = p23r.clone
+          o23cmz = p23c.clone
+          o23rmz = p23r.clone
+        o23cpz += p23z
+        o23rpz += p23z
+        o23cmz -= p23z
+        o23rmz -= p23z
+        check(o23cpz.asVector == v23cpz)
+        check(o23rpz.asVector == v23rpz)
+        check(o23cmz.asVector == v23cmz)
+        check(o23rmz.asVector == v23rmz)
+
+      test "order and plus minus of a complex matrix (float32)":
+        let
+          p23c32 = p23c.to32
+          p23r32 = p23r.to32
+          v23m32 = v23m.to32
+          v23c32 = v23c.to32
+          v23r32 = v23r.to32
+          v23v32 = v23v.to32
+
+          p23z32 = p23z.to32
+          v23cpz32 = v23cpz.to32
+          v23rpz32 = v23rpz.to32
+          v23cmz32 = v23cmz.to32
+          v23rmz32 = v23rmz.to32
+
+        check(p23c32.asVector == v23m32)
+        check(p23r32.asVector == v23m32)
+        check(not ((p23c32 + p23r32).asVector == v23c32))
+        check((p23c32 + p23r32).asVector == v23r32)
+        check((p23r32 + p23c32).asVector == v23v32)
+        check((p23r32 + p23r32).asVector == v23c32)
+        check((p23c32 + p23c32).asVector == v23c32)
+
+        var
+          o23cr32 = p23c32.clone
+          o23rc32 = p23r32.clone
+          o23rr32 = p23r32.clone
+          o23cc32 = p23c32.clone
+        o23cr32 += p23r32
+        o23rc32 += p23c32
+        o23rr32 += p23r32
+        o23cc32 += p23c32
+        check(o23cr32.asVector == v23r32)
+        check(o23rc32.asVector == v23v32)
+        check(o23rr32.asVector == v23c32)
+        check(o23cc32.asVector == v23c32)
+
+        check((p23c32 + p23z32).asVector == v23cpz32)
+        check((p23r32 + p23z32).asVector == v23rpz32)
+        check((p23c32 - p23z32).asVector == v23cmz32)
+        check((p23r32 - p23z32).asVector == v23rmz32)
+
+        var
+          o23cpz32 = p23c32.clone
+          o23rpz32 = p23r32.clone
+          o23cmz32 = p23c32.clone
+          o23rmz32 = p23r32.clone
+        o23cpz32 += p23z32
+        o23rpz32 += p23z32
+        o23cmz32 -= p23z32
+        o23rmz32 -= p23z32
+        check(o23cpz32.asVector == v23cpz32)
+        check(o23rpz32.asVector == v23rpz32)
+        check(o23cmz32.asVector == v23cmz32)
+        check(o23rmz32.asVector == v23rmz32)
+
+run()
diff --git a/tests/cla/cdet.nim b/tests/cla/cdet.nim
new file mode 100644
index 0000000..ff61e91
--- /dev/null
+++ b/tests/cla/cdet.nim
@@ -0,0 +1,32 @@
+# complex linear algebra solvers
+
+import unittest, complex, neo/clasolve, neo/cla, neo/dense
+
+proc run() =
+  suite "determinant complex computations":
+    test "determinant of a complex matrix":
+      let
+        a = matrix(@[
+          @[1.0, 4.0],
+          @[2.0, 3.0]
+        ]).toComplex
+        b = matrix(@[
+          @[0.0, -1.0],
+          @[1.0, 0.0]
+        ]).toComplex
+        m = matrix(@[
+          @[1.0, 0.0, 1.0],
+          @[0.0, 1.0, 2.0]
+        ]).toComplex
+
+      test "determinant of a complex matrix (float64)":
+        check((det(a) + complex(5.0, 0.0)) =~ im(0.0)) # 5 * -1
+        check((det(b) + complex(-1.0, 0.0)) =~ im(0.0)) # i * -i
+        check(det(m) =~ im(0.0)) # TODO: must be assertion error ?
+
+      test "determinant of a complex matrix (float32)":
+        check((a.to32.det + complex(5.0'f32, 0.0'f32)) =~ im(0.0'f32)) # 5 * -1
+        check((b.to32.det + complex(-1.0'f32, 0.0'f32)) =~ im(0.0'f32)) # i * -i
+        check(m.to32.det =~ im(0.0'f32)) # TODO: must be assertion error ?
+
+run()
diff --git a/tests/cla/ceig.nim b/tests/cla/ceig.nim
new file mode 100644
index 0000000..21d1d56
--- /dev/null
+++ b/tests/cla/ceig.nim
@@ -0,0 +1,139 @@
+# complex linear algebra solvers
+
+import unittest, complex, neo/clasolve, neo/cla, neo/dense
+import math
+import strformat, strutils
+
+proc fmtM*[T](a: Matrix[Complex[T]]): string=
+  var r = @["[\n"]
+  for i in 0..<a.M:
+    r.add("[")
+    for j in 0..<a.N: r.add(fmt" ({a[i, j].re:>21.17f}, {a[i, j].im:>21.17f})")
+    r.add("]\n")
+  r.add("]")
+  result = r.join
+
+proc fmtM*[T: SomeFloat](a: Matrix[T]): string=
+  var r = @["[\n"]
+  for i in 0..<a.M:
+    r.add("[")
+    for j in 0..<a.N: r.add(fmt" {a[i, j]:>21.17f}")
+    r.add("]\n")
+  r.add("]")
+  result = r.join
+
+proc checkEvsEvecs*[T](a: Matrix[Complex[T]],
+  aev: seq[complex.Complex[T]], aevec: seq[Vector[Complex[T]]]): bool=
+  let (evs, evecs) = a.geneig
+  for j in 0..<a.N:
+    check(evs.real[j].toComplex(evs.img[j]) =~ aev[j])
+    check(evecs.real[j].toComplex(evecs.img[j]) =~ aevec[j])
+  result = true
+
+proc checkLambdas*[T](a: Matrix[Complex[T]],
+  aev: seq[complex.Complex[T]], aevec: seq[Vector[Complex[T]]]): bool=
+  let (_, S, Vh) = a.usvd # (U, S, Vh)
+  # TODO: verification by some calculations
+  when false:
+    block:
+      let (U, S, Vh) = a.usvd
+      echo fmt"A = {a.fmtM}"
+      echo fmt"U = {U.fmtM}"
+      echo fmt"S = {S.fmtM}"
+      echo fmt"Vh = {Vh.fmtM}"
+      echo fmt"Vh * S.vdiag.diag * Vh.inv = {(Vh * S.vdiag.diag * Vh.inv).fmtM}"
+  check(Vh * S.vdiag.diag * Vh.inv =~ a)
+  check(Vh.inv * a * Vh =~ S.vdiag.diag)
+  result = true
+
+proc run() =
+  suite "eigenvalues eigenvectors complex computations":
+    const
+      r1div5 = sqrt(5.0) / 5.0
+      r1div2 = sqrt(2.0) / 2.0
+
+    let
+      are = @[
+        @[
+          @[1.0, 4.0],
+          @[2.0, 3.0]
+        ],
+        @[ # complex eigenvalue
+          @[0.0, -1.0],
+          @[1.0, 0.0]
+        ],
+        @[ # diagonal matrix
+          @[1.0, 0.0],
+          @[0.0, 2.0]
+        ],
+        @[ # multiple eigenvalue solutions
+          @[1.0, 0.0],
+          @[0.0, 1.0]
+        ]
+      ]
+
+      aevs = @[
+        EigenValues[float64](
+          real: @[-1.0, 5.0],
+          img: @[0.0, 0.0]
+        ),
+        EigenValues[float64]( # complex eigenvalue
+          real: @[0.0, 0.0],
+          img: @[1.0, -1.0]
+        ),
+        EigenValues[float64]( # diagonal matrix
+          real: @[1.0, 2.0],
+          img: @[0.0, 0.0]
+        ),
+        EigenValues[float64]( # multiple eigenvalue solutions
+          real: @[1.0, 1.0],
+          img: @[0.0, 0.0]
+        )
+      ]
+
+      aevecs = @[
+        EigenVectors[float64](
+          real: @[vector(@[2.0 * r1div5, -r1div5]), vector(@[r1div2, r1div2])],
+          img: @[vector(@[0.0, 0.0]), vector(@[0.0, 0.0])]
+        ),
+        EigenVectors[float64]( # complex eigenvalue
+          real: @[vector(@[r1div2, 0.0]), vector(@[r1div2, 0.0])],
+          img: @[vector(@[0.0, -r1div2]), vector(@[0.0, r1div2])]
+        ),
+        EigenVectors[float64]( # diagonal matrix
+          real: @[vector(@[1.0, 0.0]), vector(@[0.0, 1.0])],
+          img: @[vector(@[0.0, 0.0]), vector(@[0.0, 0.0])]
+        ),
+        EigenVectors[float64]( # multiple eigenvalue solutions
+          real: @[vector(@[1.0, 0.0]), vector(@[0.0, 1.0])],
+          img: @[vector(@[0.0, 0.0]), vector(@[0.0, 0.0])]
+        )
+      ]
+
+    test "eigenvalues eigenvectors of a complex matrix (float64)":
+      for i in 0..<are.len:
+        let a = matrix(are[i]).toComplex
+        var aev = newSeq[complex.Complex[float64]](a.N)
+        var aevec = newSeq[Vector[Complex[float64]]](a.N)
+        for j in 0..<a.N:
+          aev[j] = aevs[i].real[j].toComplex(aevs[i].img[j])
+          aevec[j] = aevecs[i].real[j].toComplex(aevecs[i].img[j])
+        test fmt"checkEvsEvecs(f64) {i}":
+          check(checkEvsEvecs(a, aev, aevec))
+        test fmt"checkLambdas(f64) {i}":
+          check(checkLambdas(a, aev, aevec))
+
+    test "eigenvalues eigenvectors of a complex matrix (float32)":
+      for i in 0..<are.len:
+        let a32 = matrix(are[i]).toComplex.to32
+        var aev32 = newSeq[complex.Complex[float32]](a32.N)
+        var aevec32 = newSeq[Vector[Complex[float32]]](a32.N)
+        for j in 0..<a32.N:
+          aev32[j] = aevs[i].real[j].toComplex(aevs[i].img[j]).to32
+          aevec32[j] = aevecs[i].real[j].toComplex(aevecs[i].img[j]).to32
+        test fmt"checkEvsEvecs(f32) {i}":
+          check(checkEvsEvecs(a32, aev32, aevec32))
+        test fmt"checkLambdas(f32) {i}":
+          check(checkLambdas(a32, aev32, aevec32))
+
+run()
diff --git a/tests/cla/chermitian.nim b/tests/cla/chermitian.nim
new file mode 100644
index 0000000..7e53454
--- /dev/null
+++ b/tests/cla/chermitian.nim
@@ -0,0 +1,36 @@
+# complex linear algebra solvers
+
+import unittest, complex, neo/cla, neo/dense
+
+proc run() =
+  suite "hermitian matrix complex computations":
+    test "hermitian matrix":
+      let
+        cm22h = matrix(@[
+          @[complex(1.0, 0.0), im(1.0)],
+          @[im(-1.0), complex(1.0, 0.0)]]) # [[1, i], [-i, 1]]
+        cm22t = matrix(@[
+          @[complex(1.0, 0.0), im(-1.0)],
+          @[im(1.0), complex(1.0, 0.0)]]) # [[1, -i], [i, 1]]
+
+      test "hermitian matrix (float64)":
+        check(not (cm22h.t == cm22h))
+        check(not (cm22h.t =~ cm22h))
+        check(cm22h.H == cm22h)
+        check(cm22h.H =~ cm22h)
+        check(cm22t.H == cm22t)
+        check(cm22t.H =~ cm22t)
+
+      test "hermitian matrix (float32)":
+        let
+          cm22h32 = cm22h.to32
+          cm22t32 = cm22t.to32
+
+        check(not (cm22h32.t == cm22h32))
+        check(not (cm22h32.t =~ cm22h32))
+        check(cm22h32.H == cm22h32)
+        check(cm22h32.H =~ cm22h32)
+        check(cm22t32.H == cm22t32)
+        check(cm22t32.H =~ cm22t32)
+
+run()
diff --git a/tests/cla/cinv.nim b/tests/cla/cinv.nim
new file mode 100644
index 0000000..8d93854
--- /dev/null
+++ b/tests/cla/cinv.nim
@@ -0,0 +1,76 @@
+# complex linear algebra solvers
+
+import unittest, complex, neo/clasolve, neo/cla, neo/dense
+
+proc run() =
+  suite "invert complex computations":
+    test "invert of a complex matrix":
+      let
+        e = matrix(@[
+          @[1.0, 0.0],
+          @[0.0, 1.0]
+        ]).toComplex
+        a = matrix(@[
+          @[1.0, 4.0],
+          @[2.0, 3.0]
+        ]).toComplex
+        ainv = matrix(@[
+          @[-0.6, 0.8],
+          @[0.4, -0.2]
+        ]).toComplex
+        b = matrix(@[
+          @[0.0, -1.0],
+          @[1.0, 0.0]
+        ]).toComplex # same as qj
+        binv = matrix(@[
+          @[0.0, 1.0],
+          @[-1.0, 0.0]
+        ]).toComplex # same as qj.H
+        qi = matrix(@[
+          @[0.0, 1.0],
+          @[1.0, 0.0]
+        ]) * im(-1.0) # Unitary Matrix
+        qj = matrix(@[
+          @[im(0.0), im(-1.0)],
+          @[im(1.0), im(0.0)]
+        ]) * im(-1.0) # Unitary Matrix
+        qk = matrix(@[
+          @[1.0, 0.0],
+          @[0.0, -1.0]
+        ]) * im(-1.0) # Unitary Matrix
+
+      test "invert of a complex matrix (float64)":
+        check(inv(e) =~ e)
+        check(inv(a) =~ ainv)
+        check(inv(b) =~ binv)
+
+        check(qi.inv =~ qi.H)
+        check(qi * qi.H =~ e)
+
+        check(qj.inv =~ qj.H)
+        check(qj * qj.H =~ e)
+
+        check(qk.inv =~ qk.H)
+        check(qk * qk.H =~ e)
+
+      test "invert of a complex matrix (float32)":
+        let
+          e32 = e.to32
+          qi32 = qi.to32
+          qj32 = qj.to32
+          qk32 = qk.to32
+
+        check(e32.inv =~ e32)
+        check(a.to32.inv =~ ainv.to32)
+        check(b.to32.inv =~ binv.to32)
+
+        check(qi32.inv =~ qi32.H)
+        check(qi32 * qi32.H =~ e32)
+
+        check(qj32.inv =~ qj32.H)
+        check(qj32 * qj32.H =~ e32)
+
+        check(qk32.inv =~ qk32.H)
+        check(qk32 * qk32.H =~ e32)
+
+run()
diff --git a/tests/cla/cpauli.nim b/tests/cla/cpauli.nim
new file mode 100644
index 0000000..620e143
--- /dev/null
+++ b/tests/cla/cpauli.nim
@@ -0,0 +1,76 @@
+# complex linear algebra solvers
+
+import unittest, complex, neo/cla, neo/dense
+
+proc run() =
+  suite "pauli matrix complex computations":
+    test "pauli matrix":
+      const
+        n = 2
+        cz = im(0.0)
+        ci = im(1.0)
+        cmi = im(-1.0)
+        co = complex(1.0, 0.0)
+        cmo = complex(-1.0, 0.0)
+      let
+        cMatI = eye(n).toComplex
+        cMatmI = matrix(@[@[cmo, cz], @[cz, cmo]])
+
+      let
+        m22zr = zeros(n, n) # zeros real matrix
+        si = matrix(@[@[0.0, 1.0], @[1.0, 0.0]]).toComplex(m22zr)
+        sj = m22zr.toComplex(matrix(@[@[0.0, -1.0], @[1.0, 0.0]]))
+        sk = matrix(@[@[1.0, 0.0], @[0.0, -1.0]]).toComplex(m22zr)
+        qi = complex(0.0, -1.0) * si
+        qj = complex(0.0, -1.0) * sj
+        qk = complex(0.0, -1.0) * sk
+
+      test "pauli matrix (float64)":
+        check(cMatmI == cMatI * cmo)
+
+        check(si == matrix(@[@[cz, co], @[co, cz]]))
+        check(sj == matrix(@[@[cz, cmi], @[ci, cz]]))
+        check(sk == matrix(@[@[co, cz], @[cz, cmo]]))
+
+        check(qi == matrix(@[@[cz, cmi], @[cmi, cz]])) # == qI
+        check(qj == matrix(@[@[cz, cmo], @[co, cz]])) # == qJ
+        check(qk == matrix(@[@[cmi, cz], @[cz, ci]])) # == qK
+        check(qi * qi == cMatmI) # == -I
+        check(qj * qj == cMatmI) # == -I
+        check(qk * qk == cMatmI) # == -I
+        check(qi * qj == qk)
+        check(qj * qk == qi)
+        check(qk * qi == qj)
+        check(qj * qi == qk * cmo)
+        check(qk * qj == qi * cmo)
+        check(qi * qk == qj * cmo)
+
+      test "pauli matrix (float32)":
+        check(cMatmI.to32 == (cMatI * cmo).to32)
+
+        let
+          si32 = si.to32
+          sj32 = sj.to32
+          sk32 = sk.to32
+          qi32 = complex(0.0'f32, -1.0'f32) * si32
+          qj32 = complex(0.0'f32, -1.0'f32) * sj32
+          qk32 = complex(0.0'f32, -1.0'f32) * sk32
+
+        check(si32 == matrix(@[@[cz, co], @[co, cz]]).to32)
+        check(sj32 == matrix(@[@[cz, cmi], @[ci, cz]]).to32)
+        check(sk32 == matrix(@[@[co, cz], @[cz, cmo]]).to32)
+
+        check(qi32 == matrix(@[@[cz, cmi], @[cmi, cz]]).to32) # == qI
+        check(qj32 == matrix(@[@[cz, cmo], @[co, cz]]).to32) # == qJ
+        check(qk32 == matrix(@[@[cmi, cz], @[cz, ci]]).to32) # == qK
+        check(qi32 * qi32 == cMatmI.to32) # == -I
+        check(qj32 * qj32 == cMatmI.to32) # == -I
+        check(qk32 * qk32 == cMatmI.to32) # == -I
+        check(qi32 * qj32 == qk.to32)
+        check(qj32 * qk32 == qi.to32)
+        check(qk32 * qi32 == qj.to32)
+        check(qj32 * qi32 == (qk * cmo).to32)
+        check(qk32 * qj32 == (qi * cmo).to32)
+        check(qi32 * qk32 == (qj * cmo).to32)
+
+run()
diff --git a/tests/cla/csolve.nim b/tests/cla/csolve.nim
new file mode 100644
index 0000000..ba1ee71
--- /dev/null
+++ b/tests/cla/csolve.nim
@@ -0,0 +1,55 @@
+# complex linear algebra solvers
+
+import unittest, complex, neo/clasolve, neo/cla, neo/dense
+
+proc run() =
+  suite "solve complex computations":
+    let
+      a = matrix(@[
+        @[1.0, 4.0],
+        @[2.0, 3.0]
+      ]).toComplex
+      b = matrix(@[
+        @[0.0, -1.0],
+        @[1.0, 0.0]
+      ]).toComplex
+
+    test "solve of a complex matrix":
+      let
+        ma = eye(a.N) * complex(-1.0, 0.0) # (eye(a.N) * -1.0).toComplex
+        sma = matrix(@[
+          @[0.6, -0.8],
+          @[-0.4, 0.2]
+        ]).toComplex
+        mb = (eye(b.N) * -1.0).toComplex
+        smb = matrix(@[
+          @[0.0, -1.0],
+          @[1.0, 0.0]
+        ]).toComplex
+
+      test "solve of a complex matrix (float64)":
+        check(ma =~ mb)
+        check(a.solve(ma) =~ sma)
+        check(b.solve(mb) =~ smb)
+
+      test "solve of a complex matrix (float32)":
+        check(ma.to32 =~ mb.to32)
+        check(a.to32.solve(ma.to32) =~ sma.to32)
+        check(b.to32.solve(mb.to32) =~ smb.to32)
+
+    test "solve of a complex vector":
+      let
+        va = (ones(a.N) * -5.0).toComplex
+        sva = vector(@[-1.0, -1.0]).toComplex
+        vb = (ones(b.N) * -5.0).toComplex
+        svb = vector(@[-5.0, 5.0]).toComplex
+
+      test "solve of a complex vector (float64)":
+        check(a.solve(va) =~ sva)
+        check(b.solve(vb) =~ svb)
+
+      test "solve of a complex vector (float32)":
+        check(a.to32.solve(va.to32) =~ sva.to32)
+        check(b.to32.solve(vb.to32) =~ svb.to32)
+
+run()
diff --git a/tests/tcla.nim b/tests/tcla.nim
new file mode 100644
index 0000000..3023a7d
--- /dev/null
+++ b/tests/tcla.nim
@@ -0,0 +1,9 @@
+# complex linear algebra solvers
+
+{.push warning[ProveInit]: off .}
+
+import cla/cbase, cla/cpauli, cla/chermitian,
+  cla/cdet, cla/cinv, cla/csolve,
+  cla/ceig
+
+{. pop .}