diff --git a/Lilly1966/.ipynb_checkpoints/Untitled2-checkpoint.ipynb b/Lilly1966/.ipynb_checkpoints/Untitled2-checkpoint.ipynb new file mode 100644 index 0000000..bca5c3e --- /dev/null +++ b/Lilly1966/.ipynb_checkpoints/Untitled2-checkpoint.ipynb @@ -0,0 +1,498 @@ +{ + "cells": [ + { + "cell_type": "markdown", + "id": "d7f053a4", + "metadata": {}, + "source": [ + "# Dongarra et al. (1996) D2 method\n", + "$n$次チェビシェフ多項式$T_n$の$0\\leq n\\leq N-4$について, \n", + "\\begin{equation}\n", + "\\begin{split}\n", + "&\\dfrac{1}{24} \\sum^{N}_{\\substack{p=n+4\\\\ p\\equiv n\\; (\\mathrm{mod} \\; 2)}}{p\\left[p^2(p^2-4)^2-3n^2p^4+3n^4p^2-n^2\\left(n^2-4 \\right) ^2-48\\alpha ^2(p^2-n^2) \\right] a_p} -8\\alpha ^2(n+2)(n+1)a_{n+2} +\\alpha ^4c_na_n \\\\\n", + "&\\quad -\\dfrac{i\\alpha R}{2} \\left[\\sum^{N}_{p=2}{a_p\\sum^{}_{\\substack{m\\equiv p\\; (\\mathrm{mod} \\; 2) \\\\ |m|\\leq p-2 \\\\ |n-m|\\leq N}}{p(p^2-m^2)c_{|n-m|}b_{|n-m|} }} -\\alpha ^2\\sum^{}_{\\substack{|p|\\leq N \\\\ |n-p|\\leq N}}{c_{|p|}a_{|p|}c_{|n-p|}b_{|n-p|}} -\\sum^{}_{\\substack{|p|\\leq N \\\\ |n-p|\\leq N}}{c_{|p|}a_{|p|}\\sum^{N}_{\\substack{m=|n-p|+2\\\\ m+n\\equiv p\\; (\\mathrm{mod} \\; 2)}}{m\\left[m^2-(n-p)^2 \\right] b_m}} \\right] \\\\\n", + "&\\quad +\\lambda \\left\\{i\\alpha R\\sum^{N}_{\\substack{p=n+2\\\\ p\\equiv n\\; (\\mathrm{mod} \\; 2)}}{p(p^2-n^2)a_p} -i\\alpha ^3Rc_na_n \\right\\} =0.\n", + "\\end{split}\n", + " \\label{eq:Orszag1971-17-A44}\n", + "\\end{equation}\n", + "および境界条件から得られる係数の関係式:\n", + "$$\n", + "\\sum^{N_{\\mathrm{e}}}_{\\substack{p=0\\\\ p=0\\; (\\mathrm{mod}\\; 2)}}{a_{\\upsilon ,p}} =0,\\quad \n", + "\\sum^{N_{\\mathrm{e}}}_{\\substack{p=0\\\\ p=0\\; (\\mathrm{mod}\\; 2)}}{a_{\\chi ,p}} =0,\\quad \n", + "\\sum^{N_{\\mathrm{o}}}_{\\substack{p=1\\\\ p=1\\; (\\mathrm{mod}\\; 2)}}{a_{\\upsilon ,p}} =0,\\quad \n", + "\\sum^{N_{\\mathrm{o}}}_{\\substack{p=1\\\\ p=1\\; (\\mathrm{mod}\\; 2)}}{a_{\\chi ,p}} =0.\n", + "$$\n", + "求める固有値方程式は\n", + "$$\n", + "(A-B\\lambda )\\textbf{x} =\\textbf{0} \\quad \\rightarrow \\quad A\\textbf{x} =\\lambda B\\textbf{x} ,\n", + "$$\n", + "$$\n", + "\\textbf{x} \\equiv \\left[a_{\\upsilon ,0},a_{\\upsilon ,1},\\cdots ,a_{\\upsilon ,N},\\cdots a_{\\chi ,0},a_{\\chi ,1},\\cdots ,a_{\\chi ,N} \\right] ^T.\n", + "$$" + ] + }, + { + "cell_type": "code", + "execution_count": 1, + "id": "7e7d6392", + "metadata": {}, + "outputs": [ + { + "name": "stdout", + "output_type": "stream", + "text": [ + "Nemax = 56, Nomax = 55\n", + "Schur(Eig): NaN + NaN*im\n", + "Schur(Eig): 0.9999999999960233 - 297.9289987076738im\n", + "Schur(Eig): 0.999999999999985 - 16.075027523504648im\n", + "Schur(Eig): 1.000000000000001 - 4.968910627205674im\n", + "Schur(Eig): 0.9999999999999996 - 2.4121138123230654im\n", + "Schur(Eig): 1.0000000000000033 - 1.4332504539356725im\n", + "Schur(Eig): 1.0000000000000022 - 0.9593882638245351im\n", + "Schur(Eig): 1.0000000000000027 - 0.6949321310788117im\n", + "Schur(Eig): 1.000000000000002 - 0.5333033206133881im\n", + "Schur(Eig): 1.0000000000000013 - 0.4281176298572095im\n", + "Schur(Eig): 0.9999999999999991 - 0.35678772751771776im\n", + "Schur(Eig): 1.0000000000000022 - 0.30681720160791687im\n", + "Schur(Eig): 1.0000000000000004 - 0.2693151105778039im\n", + "Schur(Eig): 1.0000000000000022 - 0.23724029705851687im\n", + "Schur(Eig): 0.9999999999999987 - 0.20760882901759264im\n", + "Schur(Eig): 0.9999999999999998 - 0.17997354250384043im\n", + "Schur(Eig): 0.9999999999999999 - 0.15431256878733762im\n", + "Schur(Eig): 1.0000000000000022 - 0.13062551820406373im\n", + "Schur(Eig): 0.9999999999999993 - 0.10891238852201467im\n", + "Schur(Eig): 1.000000000000001 - 0.08917317971983175im\n", + "Schur(Eig): 1.0000000000000027 - 0.07140789179787063im\n", + "Schur(Eig): 1.0000000000000002 - 0.055616524756129104im\n", + "Schur(Eig): 1.0000000000000002 - 0.00034674011002908437im\n", + "Schur(Eig): 1.0000000000000002 - 0.0023206609902555356im\n", + "Schur(Eig): 1.0000000000000018 - 0.006268502750680677im\n", + "Schur(Eig): 0.9999999999999999 - 0.012190265391334873im\n", + "Schur(Eig): 0.9999999999999999 - 0.04179907859460324im\n", + "Schur(Eig): 1.0000000000000004 - 0.020085948912205928im\n", + "Schur(Eig): 0.9999999999999998 - 0.02995555331329627im\n", + "Schur(Eig): NaN + NaN*im\n", + "Schur(Eig): NaN + NaN*im\n", + "Schur(Eig): NaN + NaN*im\n", + "Schur(Eig): NaN + NaN*im\n", + "Schur(Eig): NaN + NaN*im\n", + "Schur(Eig): NaN + NaN*im\n", + "Schur(Eig): NaN + NaN*im\n", + "Schur(Eig): NaN + NaN*im\n", + "Schur(Eig): NaN + NaN*im\n", + "Schur(Eig): NaN + NaN*im\n", + "Schur(Eig): NaN + NaN*im\n", + "Schur(Eig): NaN + NaN*im\n", + "Schur(Eig): NaN + NaN*im\n", + "Schur(Eig): NaN + NaN*im\n", + "Schur(Eig): NaN + NaN*im\n", + "Schur(Eig): NaN + NaN*im\n", + "Schur(Eig): NaN + NaN*im\n", + "Schur(Eig): NaN + NaN*im\n", + "Schur(Eig): NaN + NaN*im\n", + "Schur(Eig): NaN + NaN*im\n", + "Schur(Eig): NaN + NaN*im\n", + "Schur(Eig): NaN + NaN*im\n", + "Schur(Eig): NaN + NaN*im\n", + "Schur(Eig): NaN + NaN*im\n", + "Schur(Eig): NaN + NaN*im\n", + "Schur(Eig): NaN + NaN*im\n", + "Schur(Eig): NaN + NaN*im\n", + "Schur(Eig): NaN + NaN*im\n" + ] + } + ], + "source": [ + "#using .Lilly1966_solver\n", + "using OffsetArrays # 配列要素番号を任意開始にするモジュール (行列計算時は 1 始まりの通常 matrix に copy する)\n", + "using LinearAlgebra\n", + "\n", + "#-- Drawing parameters \n", + "N = 56 # wavenumber in Chebyshev\n", + "Ncal = N + 1\n", + "Ntot = 2 * (N + 1)\n", + "\n", + "# initialize working arrays\n", + "# A: non-lambda coefficients, B: lambda coefficients\n", + "A_J = reshape(zeros(Complex{Float64},Ntot,Ntot),Ntot,Ntot) # Matrix for calculation\n", + "B_J = reshape(zeros(Complex{Float64},Ntot,Ntot),Ntot,Ntot) # Matrix for calculation\n", + "### --- From here, supposing matrix element as row, col as in mathematical notation\n", + "Acal = reshape(zeros(Complex{Float64},Ntot,Ntot),Ntot,Ntot) # Matrix for calculation\n", + "Bcal = reshape(zeros(Complex{Float64},Ntot,Ntot),Ntot,Ntot) # Matrix for calculation\n", + "A_vv = OffsetArray(Acal[1:Ncal,1:Ncal],0:N,0:N) # coefficient matrix for v in v Chebyshev\n", + "A_vchi = OffsetArray(Acal[1:Ncal,1:Ncal],0:N,0:N) # coefficient matrix for chi in v Chebyshev\n", + "A_chiv = OffsetArray(Acal[1:Ncal,1:Ncal],0:N,0:N) # coefficient matrix for chi in chi Chebyshev\n", + "A_chichi = OffsetArray(Acal[1:Ncal,1:Ncal],0:N,0:N) # coefficient matrix for chi in chi Chebyshev\n", + "B_chichi = OffsetArray(Acal[1:Ncal,1:Ncal],0:N,0:N) # coefficient matrix for chi in chi Chebyshev\n", + "### --- From here, supposing matrix element as col, row as in julia notation (only using these arrays for matrix calc)\n", + "b = OffsetArray(zeros(Float64,Ncal),0:N) # A vector composed of Chebyshev coefficients for the basic state flow\n", + "c = OffsetArray(zeros(Float64,Ncal),0:N) # factorized vector (c_0=2, c_n=1, (n>0))\n", + "\n", + "# setting parameters\n", + "α = 1.0\n", + "Re = 10000.0\n", + "\n", + "# setting the Chebyshev coefficients for the basic flow\n", + "b[0] = 0.5\n", + "#b[2] = -b[0]\n", + "\n", + "# setting the factorized vector\n", + "c[0] = 2.0\n", + "c[1:N] .= 1.0\n", + "\n", + "# setting maximum even and odd numbers from N\n", + "Nemax = 0\n", + "Nomax = 1\n", + "for i in 0:2:N\n", + " Nemax = i\n", + "end\n", + "for i in 1:2:N\n", + " Nomax = i\n", + "end\n", + "println(\"Nemax = \", Nemax, \", Nomax = \", Nomax)\n", + "\n", + "# setting matrices\n", + "# 支配方程式から構築される部分 (チェビシェフ次数 0<=n<=N-2 の範囲)\n", + "# この段階では境界条件で成り立つ係数の関係は考慮しない\n", + "for n in 0:N-2\n", + " n2 = n^2\n", + " \n", + " # In v rows, a_{v,n} and a_{chi,n}\n", + " for p in n+2:2:N\n", + " p2 = p^2\n", + " A_vv[n,p] = p * (p2 - n2)\n", + " end\n", + " \n", + " A_vv[n,n] = - c[n] * α^2\n", + " A_vchi[n,n] = - c[n]\n", + " \n", + " # In chi rows, a_{v,n} and a_{chi,n}\n", + " for p in n+2:2:N\n", + " p2 = p^2\n", + " A_chichi[n,p] = p * (p2 - n2)\n", + " end\n", + " \n", + " A_chichi[n,n] = - c[n] * α^2\n", + "\n", + " for p in -N:N\n", + " abs_n_p = abs(n - p)\n", + " abs_p = abs(p)\n", + " if abs_n_p <= N\n", + " A_chichi[n,abs_p] = A_chichi[n,abs_p] + complex(0.0, - α * Re * c[abs_p] * c[abs_n_p] * b[abs_n_p])\n", + " end\n", + "\n", + " if abs_n_p+2 <= N\n", + " acoef = 0.0\n", + " for m in abs_n_p+2:2:N\n", + " m2 = m^2\n", + " acoef = acoef + m * (m2 - (n - p)^2) * b[m]\n", + " end\n", + " A_chiv[n,abs_p] = A_chiv[n,abs_p] + complex(0.0,α * Re * c[abs_p] * acoef)\n", + " end\n", + " end\n", + " \n", + " # In chi rows, lambda coefficients\n", + " B_chichi[n,n] = complex(0.0, - α * Re * c[n])\n", + "\n", + "end\n", + "\n", + "# Set boundary conditions\n", + "for p in 0:2:N # even mode\n", + " A_vv[N-1,p] = 1.0\n", + " A_chichi[N-1,p] = 1.0\n", + "end\n", + "for p in 1:2:N # even mode\n", + " A_vv[N,p] = 1.0\n", + " A_chichi[N,p] = 1.0\n", + "end\n", + "\n", + "# Gather A_v, A_chi, B_chi -> Acal, Bcal\n", + "Acal[1:Ncal,1:Ncal] .= A_vv[0:N,0:N]\n", + "Acal[1:Ncal,Ncal+1:Ntot] .= A_vchi[0:N,0:N]\n", + "Acal[Ncal+1:Ntot,1:Ncal] .= A_chiv[0:N,0:N]\n", + "Acal[Ncal+1:Ntot,Ncal+1:Ntot] .= A_chichi[0:N,0:N]\n", + "Bcal[Ncal+1:Ntot,Ncal+1:Ntot] .= B_chichi[0:N,0:N]\n", + "\n", + "# Convert Transpose type to Matrix (Array) type (ドット演算して要素だけコピーしている)\n", + "A_J .= transpose(Acal) # A_J = Transpose type\n", + "B_J .= transpose(Bcal) # B_J = Transpose type\n", + "\n", + "F = schur(A_J,B_J)\n", + "D2cal = F.α ./ F.β\n", + "for i in 1:Ncal\n", + " println(\"Schur(Eig): \", D2cal[i])\n", + "end\n", + "\n", + "#G = eigvecs(Ccal)\n", + "#H = inv(G) * Ccal * G\n", + "#ncmax = findmax(F)[1]\n", + "#Ae = OffsetArray(A[0:div(N,2),0:div(N,2)],0:div(N,2),0:div(N,2))\n", + "#Ao = OffsetArray(A[1:div(N,2),1:div(N,2)],1:div(N,2),1:div(N,2))\n", + "#Be = OffsetArray(B[0:div(N,2),0:div(N,2)],0:div(N,2),0:div(N,2))\n", + "#Bo = OffsetArray(B[1:div(N,2),1:div(N,2)],1:div(N,2),1:div(N,2))\n", + "#for n in 0:div(N,2)\n", + "# for p in 0:div(N,2)\n", + "# Ae[n,p] = A[2*n,2*p]\n", + "# Be[n,p] = B[2*n,2*p]\n", + "# end\n", + "#end\n", + "#for n in 1:div(N,2)\n", + "# for p in 1:div(N,2)\n", + "# Ao[n,p] = A[2*n-1,2*p-1]\n", + "# Bo[n,p] = B[2*n-1,2*p-1]\n", + "# end\n", + "#end\n", + "#Acal[1:div(N,2)-1,1:div(N,2)-1] .= Ae[0:div(N,2)-2,0:div(N,2)-2]\n", + "#Bcal[1:div(N,2)-1,1:div(N,2)-1] .= Be[0:div(N,2)-2,0:div(N,2)-2]\n", + "#Ccal = inv(Bcal[1:div(N,2)-1,1:div(N,2)-1]) * Acal[1:div(N,2)-1,1:div(N,2)-1]\n", + "##Ccal = inv(Bcal) * Acal\n", + "#Dcal = eigvals(Ccal)\n", + "#F = imag.(Dcal)\n", + "#for i in 1:div(N,2)-1\n", + "# println(F[i])\n", + "#end\n", + "\n", + "#Acal[1:div(N,2)-2,1:div(N,2)-2] .= Ao[1:div(N,2)-2,1:div(N,2)-2]\n", + "#Bcal[1:div(N,2)-2,1:div(N,2)-2] .= Bo[1:div(N,2)-2,1:div(N,2)-2]\n", + "#Ccal = inv(Bcal[1:div(N,2)-2,1:div(N,2)-2]) * Acal[1:div(N,2)-2,1:div(N,2)-2]\n", + "##Ccal = inv(Bcal) * Acal\n", + "#Dcal = eigvals(Ccal)\n", + "#F = imag.(Dcal)\n", + "#for i in 1:div(N,2)-2\n", + "# println(F[i])\n", + "#end\n" + ] + }, + { + "cell_type": "code", + "execution_count": 2, + "id": "1255c473", + "metadata": {}, + "outputs": [ + { + "ename": "UndefVarError", + "evalue": "\u001b[91mUndefVarError: Ncal not defined\u001b[39m", + "output_type": "error", + "traceback": [ + "\u001b[91mUndefVarError: Ncal not defined\u001b[39m", + "", + "Stacktrace:", + " [1] top-level scope at In[2]:9" + ] + } + ], + "source": [ + "#-- Drawing\n", + "##########\n", + "# Plot #\n", + "##########\n", + "using PyPlot\n", + "\n", + "rc(\"font\", family=\"IPAPGothic\")\n", + "fig = figure(\"pyplot_majorminor\",figsize=(7,5))\n", + "\n", + "x = real.(Dcal[3:Ncal-4])\n", + "y = imag.(Dcal[3:Ncal-4])\n", + "\n", + "ax = gca()\n", + "\n", + "cp = scatter(x,y)\n", + "\n", + "ax = gca()\n", + "\n", + "xlabel(\"X\")\n", + "ylabel(\"Y\")\n", + "#if i == 1 #div(ne,2)-1\n", + "# plt.colorbar(cp)\n", + "#end\n", + "grid(\"on\")\n", + "\n", + "fig.canvas.draw()\n", + "gcf() # Needed for IJulia to plot inline\n", + "\n", + "#PyPlot.title(\"Re = Re[i]\")\n", + "\n", + "fname = \"test.png\"\n", + "#savefig(\"B-init_2d.pdf\")\n", + "savefig(fname)\n", + "#psi\n" + ] + }, + { + "cell_type": "code", + "execution_count": 66, + "id": "623f8c38", + "metadata": {}, + "outputs": [ + { + "data": { + "text/plain": [ + "Figure(PyObject
)" + ] + }, + "metadata": {}, + "output_type": "display_data" + } + ], + "source": [ + "#-- Drawing\n", + "##########\n", + "# Plot #\n", + "##########\n", + "#using PyPlot\n", + "\n", + "#rc(\"font\", family=\"IPAPGothic\")\n", + "#fig = figure(\"pyplot_majorminor\",figsize=(7,5))\n", + "\n", + "nx = Ncal-4\n", + "ny = Ncal-4\n", + "\n", + "xax = reshape(zeros(nx,ny),nx,ny)\n", + "yax = reshape(zeros(nx,ny),nx,ny)\n", + "\n", + "for i in 1:nx\n", + " xax[i,1:ny] .= i # xmin + (xmax-xmin) * (i-1) / (nx-1)\n", + "end\n", + "for i in 1:ny\n", + " yax[1:nx,i] .= i # ymin + (ymax-ymin) * (i-1) / (ny-1)\n", + "# yax[1:nx,i] .= Re[i] # ymin + (ymax-ymin) * (i-1) / (ny-1)\n", + "end\n", + "\n", + "cont_val = reshape(zeros(nx,ny),nx,ny)\n", + "shade_val = reshape(zeros(nx,ny),nx,ny)\n", + "#for i in 1:nr#div(nr,2)-1:div(nr,2)+1\n", + " cont_val = imag.(Ccal[1:nx,1:ny])\n", + " #for j in 1:ne\n", + " shade_val = imag.(A_J[1:nx,1:ny])\n", + " #end\n", + "#draw_num2 = 2\n", + "#if draw_num2 == 1\n", + "# cp = contourf(xax[1:nx,1:nx], xax'[1:nx,1:nx], Pftime[1:nx,1:nx,1], levels=[-15.0, -10.0, -5.0, 0.0, 5.0, 10.0, 15.0, 20.0], origin=\"image\", cmap=ColorMap(\"viridis\"), extend=\"both\")\n", + " #cp = contourf(xax[1:nx,1:nx], xax'[1:nx,1:nx], Pftime[1:nx,1:nx,2], levels=[-2.0, -1.0, -0.5, -0.25, 0.25, 0.5, 1.0, 1.5, 2.0], origin=\"image\", cmap=ColorMap(\"viridis\"), extend=\"both\")\n", + " #cp = contourf(xax[1:nx,1:no], xax'[1:nx,1:no], Pftime[1:nx,1:nx,2]*Hop'*inv(Hop*Pftime[1:nx,1:nx,2]*Hop'+Ro), origin=\"image\", cmap=ColorMap(\"viridis\"), extend=\"both\")\n", + "#elseif draw_num2 == 2\n", + " ##cp = contour(xax[1:nx,1:ny], yax[1:nx,1:ny], shade_val[1:nx,1:ny], 8, levels=[-0.01, -0.005, -0.003, -0.002, -0.001, -0.0005, 0.0, 0.0005, 0.001, 0.002, 0.003, 0.005, 0.01])\n", + " cp = contourf(xax[1:nx,1:ny], yax[1:nx,1:ny], cont_val[1:nx,1:ny], cmap=ColorMap(\"viridis\"), extend=\"both\")\n", + " ##cp = contourf(xax[1:nx,1:ny], yax[1:nx,1:ny], cont_val[1:nx,1:ny], 8, levels=[-1.0, -0.5, -0.2, -0.1, 0.0, 0.1, 0.2, 0.5, 1.0], cmap=ColorMap(\"viridis\"), extend=\"both\")\n", + " #cp = contourf(xax[1:nx,1:ny], yax[1:nx,1:ny], shade_val[1:nx,1:ny], 10, cmap=ColorMap(\"viridis\"), extend=\"both\")\n", + "#elseif draw_num2 == 3\n", + "# cp = contourf(xoax[1:nx,1:nto-1], toax[1:nx,1:nto-1], evec_max[1:nx,1:nto-1], 10, cmap=ColorMap(\"viridis\"), extend=\"both\")\n", + "#elseif draw_num2 == 4\n", + "# cp = contourf(xoax[1:nx,1:40], toax[1:nx,1:40], x_inc[1:nx,1:40], levels=[-2.0, -1.0, -0.5, -0.25, 0.25, 0.5, 1.0, 1.5, 2.0], cmap=ColorMap(\"viridis\"), extend=\"both\")\n", + "#elseif draw_num2 == 5\n", + "# cp = contourf(xoax[1:no,1:40], toax[1:no,1:40], y_innov[1:no,1:40], levels=[-2.0, -1.0, -0.5, -0.25, 0.25, 0.5, 1.0, 1.5, 2.0], cmap=ColorMap(\"viridis\"), extend=\"both\")\n", + "#end\n", + "#ax.label(cp, inline=1, fontsize=10)\n", + "#legend()\n", + "ax = gca()\n", + "\n", + "xlabel(\"X\")\n", + "ylabel(\"Y\")\n", + "#if i == 1 #div(ne,2)-1\n", + " plt.colorbar(cp)\n", + "#end\n", + "grid(\"on\")\n", + "\n", + "PyPlot.title(\"Re = Re[i]\")\n", + "\n", + "#########################\n", + "# Set tick dimensions #\n", + "#########################\n", + "#ax.xaxis.set_tick_params(which=\"major\",length=5,width=2,labelsize=10)\n", + "#ax.xaxis.set_tick_params(which=\"minor\",length=5,width=2)\n", + "\n", + "fig.canvas.draw() # Update the figure\n", + "gcf() # Needed for IJulia to plot inline\n", + "#fname = \"Re\" * string(1) * \".png\"\n", + "fname = \"test.png\"\n", + "#savefig(\"B-init_2d.pdf\")\n", + "savefig(fname)\n", + "#psi\n", + "#end\n", + "#println(shade_val[1:na,1])" + ] + }, + { + "cell_type": "code", + "execution_count": 47, + "id": "af6a72bd", + "metadata": {}, + "outputs": [ + { + "name": "stdout", + "output_type": "stream", + "text": [ + "ComplexF64[0.0 + 0.0im 0.0 + 0.0im; 0.0 + 0.0im 0.0 + 0.0im]\n" + ] + }, + { + "data": { + "text/plain": [ + "GeneralizedSchur{ComplexF64, Matrix{ComplexF64}}\n", + "S factor:\n", + "2×2 Matrix{ComplexF64}:\n", + " 2.23607-0.0im 0.0+0.0im\n", + " 0.0+0.0im 2.23607+0.0im\n", + "T factor:\n", + "2×2 Matrix{ComplexF64}:\n", + " 2.23607+0.0im 0.0+0.0im\n", + " 0.0+0.0im 2.23607+0.0im\n", + "Q factor:\n", + "2×2 Matrix{ComplexF64}:\n", + " -0.447214-0.894427im 0.0+0.0im\n", + " 0.0+0.0im 1.0+0.0im\n", + "Z factor:\n", + "2×2 Matrix{ComplexF64}:\n", + " -1.0-0.0im 0.0+0.0im\n", + " -0.0-0.0im 0.894427+0.447214im\n", + "α:\n", + "2-element Vector{ComplexF64}:\n", + " 2.23606797749979 - 0.0im\n", + " 2.23606797749979 + 0.0im\n", + "β:\n", + "2-element Vector{ComplexF64}:\n", + " 2.23606797749979 + 0.0im\n", + " 2.23606797749979 + 0.0im" + ] + }, + "execution_count": 47, + "metadata": {}, + "output_type": "execute_result" + } + ], + "source": [ + "a33 = reshape(zeros(Complex{Float64},2,2),2,2)\n", + "c33 = reshape(zeros(Complex{Float64},2,2),2,2)\n", + "println(a33)\n", + "a33[1,1] = 1.0 + 2.0im\n", + "a33[2,2] = 2.0 - 1.0im\n", + "b33 = transpose(a33)\n", + "c33 .= b33\n", + "schur(a33,c33)" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "id": "b28e1119", + "metadata": {}, + "outputs": [], + "source": [] + } + ], + "metadata": { + "kernelspec": { + "display_name": "Julia 1.6.1", + "language": "julia", + "name": "julia-1.6" + }, + "language_info": { + "file_extension": ".jl", + "mimetype": "application/julia", + "name": "julia", + "version": "1.6.1" + } + }, + "nbformat": 4, + "nbformat_minor": 5 +}