! DOUBLE Precision.
! Eigensystem is printed in 'g24.16' format.
!
! Main - for QR
! With Householder - Reduction to upper triangular form
! Symmetric or Not_Symmetric matrices
!
! Data File (Input):
! -----------------
! - Title
! - n, kodA, kodPrintA
! If kodA <0:
!	- job: parameter for "agen";
!
! If kodA =0, 1
!   - Matrix A (on rows)
!
! [If kodA >=2: matrix read from file "Fname_2". No data input. Fname_2 is selected via the standard dialog window.]
!
! - Balance, kod_Balance
! If kod_Balance /=0:
!	- diags(i), i=1,n
! - Form
! - Eps, Eps_lambda, Lnit, kodIter
! - Shift (S,s... or other chars), kodShift
! - kodVP, kod_PrintVP
! If kodVP /=0:
! - Refine (logical)
!	- If(refine): is (integer)
! - kodNorm (if kodVP /=0)
! - Check_text (C,c,P,p, F,f; or, other chars)
! - Do_Gen_Files (logical:T,F)
!
! Data meanning:
! -------------
! Title: text data (1-80 chars). Problem description.
!
! n: matrix order
! kodA >=0: matrix read from data file; kodA  <0: generated matrix (by Agen);
!	Do Input:
!	kodA =0:  upper triangle;
!	kodA =1:  full input
!
!   kodA >=2: matrix read from file "fname_2", namely:
!		kodA =2:  upper triangle.
!		kodA =3:  full matrix
!		kodA >=4: non-zero matrix elements: i, j, a(i,j).
!	 
!	Note: "Fname_2" structure:
!         - The first line is a text-line (max. 80 chars; e.g.: problem title).
!         - Matrix entries, according to kodA.
!
! kodA<0: job: Character(20): Parameter to be transmitted to "agen".
!
! kodPrintA:
!		>0: print A (given & block-tri reduced)
!		=0: print A block-tri reduced
!		<0: print nothing
!
! Balance: text data (1-20 chars) - First character is significant:
!		"B", "b": do balance. Other chars: do not balance.
! kod_Balance: =0: Balance matrix by subroutine "Balance"; 
!			   #0: Direct Balance with "diags" (subroutine "Direct_Balance").
!              Any integer if no balance.
!	diags(1:n): Diagonal elements for direct balance ( A -> D^(-1)AD )
!       
! Form = Preliminary transformation of the matrix.
!       Text data (1-20 chars) - First character is significant:
!		- "A", "a" (Auto): Symmetric matrix: Tridiagonal; Not Symmetric: Hessenberg
!		- "T", "t" : Tridiagonal
!		- "H", "h": Hessenberg
!		- Other chars.: Given A (No preliminary transformations).
!
! Eps = Tolerance to off-diagonal elements
! Eps_lambda = Tolerance to two succesive lambdas
!
! Lnit = Limiting number of QR iterations
!
! kodIter =0/#0: do / do not: print iterations
!
! Shift: text data (1-20 chars). First character is significant:
!		- "S","s"/ Other chars: with/without shift;
! kodShift =0: AK(n,n)
! kodShift >0: Eigenvalue of AK(n-1:n, n-1:n) that is closest to AK(n,n)
! kodShift <0: Minimum diagonal AK(i,i). Best for certain symmetric matrices.
!
! kodVP = #0/0: do / do not: calculate eigenvectors
! kod_PrintVP >0: Do Write eigenvectors; kod_PrintVP <=0: Do Not Write eigenvectors
!
! Refine (logical): T/F = Refine/Not refine - eigensystem by Inverse Power.
!					Note: If(refine): eigenvectors are already normalized by the 2-norm.
!		is: Defines the shift for Inverse Power:
!			is >0: shift = lambda -one unit in the "is"th significant digit of lambda.
!			Recommendation: is =4; 5,6.
!			is <0: shift =lambda -1D-4
! kodNorm =Norm code:
!		<0: do not normalize
!		=0: maximum norm (Infinity-Norm)
!		=1: 1-norm
!		=2: euclidian norm (2-Norm)
!		=3: "eng. norm": First coordinate =1.
!
! Check_text: text data (1-20 chars). First character is significant:
!		- "C","c","P","p", "F", "f"/Other chars: with/without check on eigensystem;
!		  If "F" or "f": Full printing: Differencies (A*z -lambda*z) are printed, too.
!
! Do_Gen_Files: Generate or no (T/F): Eigenvalues File & Eigenvectors File.
!
!
! (do_balance, do_Shift, do_check, full_print_check: Internal use.)
!
! Output Data:
! -----------
! ibloc(1:n): Index of block row. Internal use. If kodVP/=0: Output - as "i_flag".  .
! ibloc(i) =0/2 : "i" is the First row/Second Row of the block (i_flag =1/2).
! ibloc(i) =1 : No block (i_flag =0).
!
! If do_Balance: Ratio = maximum of: |Row_Norm -Column_Norm|/Row_Norm
! 
! Eigenvaluea are reported as multiple, if they coincid up to single precision.
!
! Data Files:
! -----------
! - stored in "Dat" Folder
! - selected via OpenFileName
!
Program QR_Method
use work
use dflib

character($maxpath) fname, fname_2

character(80) titlu, title_2

character(20) shift, check_text, ch1*1
logical do_shift, do_check, full_print_check

integer width

character(20) form, Form1*1
logical auto, tridiag, hess, given
logical symm

logical BlockReduced, real_val, TriReduced

! In "work":
!real(8), allocatable :: VEC(:,:)
!real(8), allocatable :: y(:)

!complex(8), allocatable :: VEC_Z(:,:)
!complex(8), allocatable :: y_z(:)     ! in "work"

real(8) value; complex(8) value_z;

real(8) s, lambda1
real(8), allocatable :: as(:,:), w0(:)
logical refine

character(20) job

character(10) balance_text
logical do_balance

integer status
real(8) ulp

!*
real(8) max_sub_Sub, eps_2
!*

logical Do_Gen_Files

real(4) value_4, lambda_4, value_4i, lambda_4i
real(8) re_val, imag_val


Do_Gen_Files =.FALSE.
auto =.false.; tridiag =.false.; hess =.false.; given =.false.;

call GetFile(fname, status, 'Select Input File'C)

if(status <=0) then
write(*,*) 'Error in selecting file !'
GOTO 600
endif

! Open:
open(1, file =fname, status ='old', form ='formatted')

! READ:
read(1, '(a)') titlu
read(1,*) n, kodA, kodPrintA

! Allocations:
allocate(A(n,n), Q(n,n))		!, PK(n,n))
allocate (v(n), p(n))

allocate (Q1(n,n))

A =0.
call unit2(n, Q)		! Q=I
call unit2(n, Q1)	! Q1=I

! Read/Generate:
kodA_Select: SELECT CASE(kodA)

CASE(:-1)		! kodA <0: Generate
job(1:20) =' '
read(1,*) job
call agen(A, job)

CASE(0)			! kodA =0: Read from input file
do i =1,n
read(1,*) (a(i,j), j =i,n)	! upper triangle on rows
enddo
! Fill symmetric:
do i=2,n
do j=1, i-1
a(i,j) =a(j,i)
enddo
enddo

CASE(1:)		! kodA >0: Read from file
	if(kodA ==1) then	! Input File
	do i =1,n
	read(1,*) (a(i,j), j =1,n)	! full input on rows
	enddo

	else	! Specified File "fname_2"
	call GetFile(fname_2, status, 'Select Matrix File'C)

		if(status <=0) then
		write(*,*) 'Error in selecting file !'
		GOTO 600
		endif

	open(2, file =fname_2, status ='old', form ='formatted')
	read(2,*) title_2

		if(kodA ==2) then	! read upper triangle
		do i =1,n
		read(2,*) (a(i,j), j =i,n)	! upper triangle on rows
		enddo

		! Fill symmetric:
		do i=2,n
		do j=1, i-1
		a(i,j) =a(j,i)
		enddo
		enddo

		elseif(kodA ==3) then	! full input
		do i =1,n
		read(2,*) (a(i,j), j =1,n)	! full input on rows
		enddo

		else	! kodA.=4

		40 read(2,*, err =50, end =50) i, j, a(i,j)		! non-zero entries
		goto 40
		endif

	close(2)
	endif

END SELECT kodA_Select

50 continue
! Check for symmetric matrix:
if(kodA ==0) then
simetric =.true.

else
simetric =symm(n,a)
endif

! Balance:
read(1, *) balance_text, kod_Balance
do_balance =(balance_text(1:1) =='B' .or. balance_text(1:1) =='b')

IF(do_balance) THEN
allocate(diag_bal(n)); d =1.

if(kod_Balance /=0) then
read(1,*) (diag_bal(i), i=1,n)
endif

ENDIF

read(1, '(a)') Form
Form1 =Form(1:1)

! Establish form:
select case(Form1)
case("A", "a")
auto =.true.

case("T", "t")
tridiag =.true.

case("H", "h")
hess =.true.

case default
given =.true.

end select

IF(auto) THEN
	if(simetric) then
	tridiag =.true.
	else
	hess =.true.
	endif
ENDIF

read(1,*) eps, eps_lambda, lnit, kodIter

! Store
lnit0 =lnit

read(1,*) Shift, kodShift
do_shift =shift(1:1) =='S' .or. shift(1:1) =='s'

read(1,*) kodVP, kod_PrintVP

if(kodVP /=0) then
read(1,*) refine

if(refine) then
read(1,*) is
endif

read(1,*) kodNorm

read(1,*) check_text
ch1 =check_text(1:1)
do_check = ch1 =='C' .or. ch1 =='c' .or. ch1 =='P' .or. ch1 =='p' &
			.or. ch1 =='F' .or. ch1 =='f'
full_print_check =ch1 =='F' .or. ch1 =='f'
endif

! Eigen Files
read(1,*, err =130) Do_Gen_Files

130 continue

! For printing:
n8 =min(n,8)

! Printing:
Print_A1: IF(kodPrintA >0) THEN

if(simetric) then
write(1, '(/,a)') ' Matrix A (Symmetric)'
else
write(1, '(/,a)') ' Matrix A (Not Symmetric)'
endif

call PrintA(n,A)

ENDIF Print_A1

! Balance with D: A = D^(-1)*A*D
IF(do_balance) THEN
allocate(a0(n,n))
a0 =a

if(kod_Balance ==0) then
call Balance(n, a, diag_bal, 0)

else
call Direct_Balance(n, a, diag_bal)
endif

	if(kodPrintA >=0) then
	write(1,'(/,a)') 'A-balanced'
	call PrintA(n,A)
	endif
ENDIF


! Save A - for eigenvectors check:
IF(kodVP /=0) THEN
if(.not. do_balance) then
allocate(a0(n,n))
a0 =a
endif

ENDIF


IF(tridiag) THEN
call Tridiagonal(0)

ELSEIF(hess) THEN
call Hessenberg(0)

ELSE
! Given: do nothing
ENDIF

! Store Q in Q1
Q1 =Q		!

Print_A2: IF(kodPrintA >=0) THEN

IF(tridiag) THEN
write(1, '(/, a)') ' Matrix A - Tridiagonal form'

ELSEIF(hess) THEN
write(1, '(/, a)') ' Matrix A - Hessenberg form'

ELSE
! Given:
write(1, '(/, a)') ' Matrix A - Given form'
ENDIF

call PrintA(n,A)

ENDIF Print_A2

!*
!eps_2 =max_sub_Sub(jmax)
!write(1, '(a, 1pg15.7, a, I<iwl>,a)') '   Maximum sub_Sub-Diagonal (in modulus) is:', sngl(eps_2), ' (Column ', jmax, ').'
!*

! QR:
if(n >=50) then
s =secnds(0.0)
endif

allocate(lambda(n))
allocate(lambda_z(n))
lambda =0
lambda_z =0

! Prepare:
! Store "diag":
allocate(diag(n))
do i =1,n
diag(i) =a(i,i)
enddo


! Seek blocks: Set "ibloc":
allocate(ibloc(n))
ibloc =1

iter =lnit
call QR_Z(iter, do_shift, kodShift, kodIter, BlockReduced, real_val)

if(n >=50) then
s =secnds(s)
!write(1,'(a, f7.2, a)') ' ** QR Iteration: Running Time =', s, ' sec.'
endif

!write(1,*) 'Real Eigenvalues:', real_val
write(1, '(/,a, L)') ' Real Eigenvalues:', real_val


if(kodIter ==0) then
iw =width(iter)
write(1, '(/, a, I<iw>,a)') ' Iteration ', iter, ':'
else
write(1,*)
endif

TriReduced =All(ibloc ==1)

IF(BlockReduced) THEN
	if(TriReduced) then
	write(1,*) 'Reduced Matrix - Triangular Form'

	else
	write(1,*) 'Reduced Matrix - Block-Triangular Form'
	endif

ELSE
write(1,*) 'Reduced Matrix (Intended to be Block-Triangular)'
ENDIF

Print_A3: IF(kodPrintA >=0) THEN
do i =1,n
	if(n >=20) then
	iwi =width(i)
	write(1, '(I<iwi>)') i
	endif
write(1, '( <n8>(1pg15.7) )') (A(i,j), j =1,n)
if(ibloc(i) ==2 .or. ibloc(i) ==1) write(1,*)
enddo

IF(.not.TriReduced) THEN
iwn =width(n)+1
it =iwn+5
write(1, '(/, a, t<it>, a)') ' Row', 'Block Flag'

do l=1,n
if(ibloc(l) ==0) then
i_flag =1
elseif(ibloc(l) ==1) then
i_flag =0
else
i_flag =2
endif

write(1, '(I<iwn>, t<it>, I1)') l, i_flag	!ibloc(l)
enddo
ENDIF

ENDIF Print_A3


Not_Refine_IF: IF(.not. refine) THEN

write(1, '(/,/, a)') ' Eigenvalues'
write(*, '(/, a)') ' Eigenvalues'

IF(real_val) THEN	
! Write lambda
lambda =real(lambda_z)

do i=1,n
iw =width(i)+1

write(1, '(I<iw>, 1pG24.16)') i, lambda(i)
write(*, '(I<iw>, 1pG20.10)') i, lambda(i)

enddo

ELSE
! Write lambda_z
do i =1,n
iw =width(i)+1

if(dimag(lambda_z(i)) /=0) then

write(1, '(I<iw>, 1pg24.16, 1pg24.16)') i, lambda_z(i)
write(*, '(I<iw>, 1pg20.10, 1pg20.10)') i, lambda_z(i)

else
write(1, '(I<iw>, 1pG24.16)') i, real(lambda_z(i))
write(*, '(I<iw>, 1pG20.10)') i, real(lambda_z(i))
endif

enddo

ENDIF

ENDIF Not_Refine_IF


VEC_IF: IF(kodVP /=0) THEN
deallocate(Q)			!, PK)


REAL_VAL_IF: IF( real_val ) THEN
lambda =real(lambda_z)

! Check for equal eigenvalues:
allocate(multiple(n), done(n))
multiple =1

iReal_do: do i =1,n
if(multiple(i) >1) cycle iReal_do

value =lambda(i)
value_4 =sngl( value )

	kr_do: do k =1,n
	if(k ==i) cycle kr_do

	lambda_4 =sngl( lambda(k) )

		if(lambda_4 ==value_4) then
		multiple(i) =multiple(i) +1
		endif
	
	enddo kr_do

enddo iReal_do


max_multiple =maxval(multiple)
! Write Multiple:
Write_Multiple: IF(max_multiple >1) THEN
write(1, '(/, a)') ' Multiplicity (algebraic)'

done =0

iReal2_do: do i=1,n

if(done(i) >0) cycle iReal2_do

value =lambda(i)
value_4 =sngl( value )

m_i =multiple(i)

	if(m_i >1) then
	iwi =width(i)+1; iwmi =width(m_i);
	write(1, '(I<iwi>, 1pg24.16, a, I<iwmi>)'), i, value, ': ', m_i

	done(i) =1
	endif

	kr2_do: do k =i+1,n
	lambda_4 =sngl( lambda(k) )
	
		if(lambda_4 ==value_4) then
		done(k) =2
		endif

	enddo kr2_do

enddo iReal2_do

ENDIF Write_Multiple


allocate(t(n,n))
t =a

allocate(y(n))
y =0

! Re-store diags - for the reduced matrix:
do i =1,n
diag(i) =a(i,i)				!t(i,i)
enddo

allocate(VEC(n,n))
VEC =0.

call Vec_BlocTri(VEC, kodNorm)

VEC =matmul(Q1, VEC)

! Compute eigenvectors - if balanced:
if(do_balance) then
jb_do: do j =1,n

ib_do: do i =1,n
	if(diag_bal(i) ==1) cycle ib_do
	VEC(i,j) =VEC(i,j)*diag_bal(i)
enddo ib_do

enddo jb_do
endif

Refine_IF: IF(refine) THEN
allocate(iter_refine(n))
iter_refine =0
p =0.	! To store test_val
! Restore
lnit =lnit0

if(n >50) write(*, '(/, a)') ' * Refining ...'
write(1, '(//, a)') ' * Refinement:'

!!TOL =1D-15
!!TOL =2.220446049250313E-016
!TOL =5D-16

!write(1, '(/, a, 1pg20.10)') ' TOL =', tol

allocate(TOL(n))

do ivp =1,n
TOL(ivp) =ulp(lambda(ivp))

! Skip if lambda =0:
if(lambda(ivp) ==0) cycle

call Power_Refine(ivp)
enddo



iwlnit =width(lnit)+2; iwn =width(n)+1;
it1 =iwn+4; it2 =it1+16; it3 =it1+12; it4 =it3+23;
write(1,'(t<it1>, a, t<it2>, a, t<it4>, a)') 'Iterations', 'Test Value', ' Tolerance'

do ivp =1,n
if(iter_refine(ivp) >=0) then
write(1,'(I<iwn>, I<iwlnit>, t<it3>, 1pg20.10, t<it4>, 1pg10.3)') ivp, iter_Refine(ivp), p(ivp), TOL(ivp)
else
write(1,'(I<iwn>, I<iwlnit>, t<it3>, 1pg20.10, a)') ivp, -iter_Refine(ivp), p(ivp), '(stationary)'
endif
enddo


! Write eigenvalues:
write(1, '(/,a)') ' Eigenvalues'
write(*, '(/,a)') ' Eigenvalues'

do ivp =1,n
iwi =width(ivp)

write(1, '(1x, I<iwi>, 1pg24.16)') ivp, lambda(ivp)
write(*, '(1x, I<iwi>, 1pg20.10)') ivp, lambda(ivp)
enddo

ENDIF Refine_IF


IF(kodNorm >=0) THEN
if(.not. refine .or. (refine .and. kodNorm /=2)) then
call Normalize_r(VEC, kodNorm)
endif
ENDIF

n6 =min(n,6)

! Write eigenvectors:
write_vectors: IF(kodVP >0 .and. kod_PrintVP >0) THEN
write(1, '(/, a, \)') ' Eigenvectors'

ii_DO: DO ii =1,n,6
n1 =ii
n2 =ii+5
if(n2>n) n2 =n

iw2 =width(n2)
write(1,'(/, I11,\)') n1
write(1,'(<n6-1>(23x,I<iw2>))') (i, i =n1+1, n2)

do i =1,n
write(1, '( <n6>(g24.16) )') VEC(i, n1:n2)
enddo

ENDDO ii_DO
ENDIF write_vectors
!kodCheck =0


ELSE		! REAL_VAL_IF		! Work in Complex. With real lambdas, too.
allocate(VEC_Z(n,n))
vec_z=(0.,0.)

allocate(t_z(n,n))
t_z =a

allocate(y_z(n))
y_z =0


! Re-store diags - for the reduced matrix:
do i =1,n
diag(i) =a(i,i)				!t_z(i,i)
enddo


! Establish multiplicity:
allocate(multiple(n), done(n))
multiple =1

iComplex_do: do i =1,n
if(multiple(i) >1) cycle iComplex_do

value_z =lambda_z(i)

re_val =dreal(value_z); imag_val =dimag(value_z);
value_4 =sngl( re_val ); value_4i =sngl( imag_val );

	kc_do: do k =1,n

	if(k ==i) cycle kc_do
	 
	value_z =lambda_z(k)

	re_val =dreal(value_z); imag_val =dimag(value_z);
	lambda_4 =sngl( re_val ); lambda_4i =sngl( imag_val );

		if(lambda_4 ==value_4 .and. lambda_4i ==value_4i) then
		multiple(i) =multiple(i) +1
		endif

	enddo kc_do

enddo iComplex_do


max_multiple =maxval(multiple)
! Write Multiple:
Write_Multiple_Z: IF(max_multiple >1) THEN
write(1, '(/, a)') ' Multiplicity (algebraic)'

done =0

iComplex2_do: do i =1,n

if(done(i) >0) cycle iComplex2_do

value_z =lambda_z(i)
re_val =dreal(value_z); imag_val =dimag(value_z);
value_4 =sngl( re_val ); value_4i =sngl( imag_val );

m_i =multiple(i)

	if(m_i >1) then
	iwi =width(i)+1; iwmi =width(m_i);
		if(dimag(value_z) ==0) then
		write(1, '(I<iwi>, 1pg20.10, a, I<iwmi>)') i, real(value_z), ': ', m_i

		else
		write(1,'(I<iwi>, 1pg20.10, 1pg20.10, a, I<iwmi>)') i, value_z, ': ', m_i
		endif
	done(i) =1
	endif



	kc2_do: do k =i+1,n

	value_z =lambda_z(k)
	re_val =dreal(value_z); imag_val =dimag(value_z);
	lambda_4 =sngl( re_val ); lambda_4i =sngl( imag_val );
	
		if(lambda_4 ==value_4 .and. lambda_4i ==value_4i) then
		done(k) =2
		endif

	enddo kc2_do

enddo iComplex2_do

ENDIF Write_Multiple_Z


! Eigenvectors:
!t =a
call Vec_BlocTri_Z(VEC_Z, kodNorm)

VEC_Z =matmul(Q1, VEC_Z)

! Compute eigenvectors - if balanced:
if(do_balance) then
jbz_do: do j =1,n

ibz_do: do i =1,n
	if(diag_bal(i)==1.) cycle ibz_do
	VEC_Z(i,j) =VEC_Z(i,j)*diag_bal(i)
enddo ibz_do

enddo jbz_do
endif


Refine_Z_IF: IF(refine) THEN
allocate(iter_refine(n))
iter_refine =0
p =0.	! To store test_val
! Restore
lnit =lnit0

if(n >50) write(*, '(/, a)') ' * Refining ...'
write(1, '(//, a)') ' * Refinement:'

!TOL =2.220446049250313D-016
!TOL =5D-16

!write(1, '(/, a, 1pg20.10)') ' TOL =', tol
allocate(TOL(n))

do ivp =1,n
TOL(ivp) =ulp(abs(lambda_z(ivp)))

! Skip if lambda =0:
if(lambda_z(ivp) ==0) cycle

call Power_Refine_Z(ivp)
enddo


iwlnit =width(lnit)+2; iwn =width(n)+1;
it1 =iwn+4; it2 =it1+16; it3 =it1+12; it4 =it3+24;
write(1,'(t<it1>, a, t<it2>, a, t<it4>, a)') 'Iterations', 'Test Value', ' Tolerance'

do ivp =1,n
if(iter_refine(ivp) >=0) then
write(1,'(I<iwn>, I<iwlnit>, t<it3>, 1pg20.10, t<it4>, 1pg10.3)') ivp, iter_Refine(ivp), p(ivp), TOL(ivp)
else
write(1,'(I<iwn>, I<iwlnit>, t<it3>, 1pg20.10, a)') ivp, -iter_Refine(ivp), p(ivp), '(stationary)'
endif
enddo


! Write eigenvalues:
write(1, '(/,a)') ' Eigenvalues'
write(*, '(/,a)') ' Eigenvalues'

do ivp =1,n
iwi =width(ivp)
if(imag(lambda_z(ivp)) /=0) then
write(1, '(1x, I<iwi>, 1pg24.16, 1pg24.16)') ivp, lambda_z(ivp)
write(*, '(1x, I<iwi>, 1pg20.10, 1pg20.10)') ivp, lambda_z(ivp)

else
write(1, '(1x, I<iwi>, 1pg24.16, 1pg24.16)') ivp, real(lambda_z(ivp))
write(*, '(1x, I<iwi>, 1pg20.10, 1pg20.10)') ivp, real(lambda_z(ivp))
endif

enddo

ENDIF Refine_Z_IF


IF(kodNorm >=0) THEN
if(.not. refine .or. (refine .and. kodNorm /=2)) then
call Normalize_Z(VEC_Z, kodNorm)
endif
ENDIF

! Write eigenvectors:
write_vectors_Z: IF(kodVP >0  .and. kod_PrintVP >0) THEN
write(1, '(/, a, \)') ' Eigenvectors'

! Real lambdas: Write groups of six vectors
	Real_If: IF(real_val) THEN

ii2_DO: DO ii =1,n,6
n1 =ii
n2 =ii+5
if(n2>n) n2 =n

iw2 =width(n2)
write(1,'(/, I11,\)') n1
write(1,'(<n6-1>(23x,I<iw2>))') (i, i =n1+1, n2)

do i =1,n
write(1, '( <n6>(g24.16) )') VEC(i,n1:n2)
enddo

ENDDO ii2_DO

	ELSE	! Complex: write vectors one by one
wr_ivp_do: do ivp =1,n
write(1, '(/, I8)') ivp

do l=1,n
	if(imag(vec_z(l, ivp)) /=0) then
	write(1,'(g24.16, g24.16)') vec_z(l, ivp)

	else
	write(1,'(g24.16)') real(vec_z(l,ivp))
	endif
enddo

enddo wr_ivp_do

	ENDIF Real_If
ENDIF write_vectors_Z

!kodCheck =1

ENDIF REAL_VAL_IF

ENDIF VEC_IF


! Generate eigen_values file & eigen_vectors file:
Gen_Eigen_Files_If: IF(Do_Gen_Files) THEN

call Gen_Eigen_Files(fname, real_val)
ENDIF Gen_Eigen_Files_If



! CHECK:
CHECK_IF: IF(do_check) THEN

call Check(real_val, full_print_check)
ENDIF CHECK_IF


500 continue
if(n >=50) then
!s =secnds(s)
write(1,'(//, a, f7.2, a)') ' ** QR Iteration: Running Time =', s, ' sec.'
endif


if(allocated(v)) deallocate(v)		!, p)
if(allocated(p)) deallocate(p)

if(allocated(A)) deallocate(A)
if(allocated(a0)) deallocate(a0)

if(allocated(t)) deallocate(t)
if(allocated(t_z)) deallocate(t_z)

if(allocated(Q)) deallocate(Q)
if(allocated(Q1)) deallocate(Q1)

if(allocated(y)) deallocate(y)
if(allocated(y_z)) deallocate(y_z)

if(allocated(VEC)) deallocate(VEC)
if(allocated(VEC_Z)) deallocate(VEC_Z)

if(allocated(iter_refine)) deallocate(iter_refine)

600 continue
end Program QR_Method



