from DDMPY import *
import os
import sys

mat    = ["A0", "A1", "A2", "A3", "A4"]
theta  = [1e2, 1e4, 1e6]
n      = 500
nk     = 1 
if(len(sys.argv) == 1):n = 500
if(len(sys.argv) == 2):n = int(sys.argv[1])
if(len(sys.argv) == 3):
    n  = int(sys.argv[1])
    nk = int(sys.argv[2])
rep    = "n_"+ str(n) + "_k_" + str(nk) + "_" + "EXP2_1"
for Imat in mat:
    wtheta = theta
    for Itheta in wtheta:
        if Itheta == 1e2: Stheta="1E2"
        if Itheta == 1e4: Stheta="1E4"
        if Itheta == 1e6: Stheta="1E6"
        A,DA,QTA = setmatrix4exp(which=Imat, n=n, k=nk, theta=Itheta)
        xexa = np.random.rand(A.shape[0],1)
        x0 = np.random.rand(A.shape[0],1)
        b    = A @ xexa 
        r0   = b - A @ x0
        T, V, filen = lanczos_RR(A, DA, QTA, nk=nk, Aname=Imat, Stheta=Stheta, L=None, b=r0, maxiter=A.shape[1], reortho=True)
        d = 0.0*T[...,0:1]
        d[0, 0] = la.norm(r0)
        y = la.solve(T, d)
        x = x0 + V @ y
        print("|---->",Imat, Itheta, nk)
        print("||x-xexa||_{A}=",np.sqrt( (x-xexa).T@A@(x-xexa) ) )
        print("||b-A@x||_{2}/||b||_{2}",la.norm(b-A@x)/la.norm(b) )
cmd = "R -q --vanilla < REXP2_1.R" + " --args " + filen
os.system(cmd)