1 
  2 % Givens rotation QR to find eigenpairs of symmetric tridiagonal in Octave
  3 % John Bryan 2025
  4 
  5 clear all;
  6 clc;
  7 format long;
  8 
  9 
 10 function A=GetA(N)
 11    rand("seed",1);
 12    a=10*rand(); #diagonal values between 0 and 10
 13    b=10*rand(); # sub and super diagonal values between 0 and 10
 14    A = diag(a*ones(1,N)) + diag(b*ones(1,N-1),1) + diag(b*ones(1,N-1),-1);
 15 endfunction
 16 
 17 
 18 function answer=CheckOffDiagaonal(a,n)
 19    # Check whether all nondiagonal entries 
 20    # absolute value is below chosen small number. 
 21    answer=0;
 22    for i=1:n
 23      i;
 24      for j=1:n
 25         j;
 26         if ((i!=j)&&(abs(a(i,j))>1e-14))
 27           answer=1;
 28         endif
 29      endfor
 30    endfor
 31 endfunction
 32 
 33 
 34 function b=sgn(p)
 35    # Octave sign function returns 0 on input 0
 36    if (p!=0)
 37       b=sign(p);
 38    else
 39       b=1;
 40    endif
 41 endfunction
 42 
 43 
 44 function [cosine,sine] = givens_rotation(a, b)
 45    # returns givens rotation cosine and sine given matrix entries a,b
 46    if b == 0;
 47       cosine = sgn(a);
 48       sine = 0;
 49    elseif (a == 0);
 50       cosine = 0;
 51       sine = -sgn(b);
 52    elseif (abs(a) > abs(b));
 53       tangent = b / a;
 54       secant = sgn(a) * sqrt(1 + tangent * tangent);
 55       cosine = 1 / secant;
 56       sine = -cosine * tangent;
 57    else
 58       cotangent = a / b;
 59       cosecant = sgn(b) * sqrt(1 + cotangent * cotangent);
 60       sine = -1 / cosecant;
 61       cosine = cotangent / cosecant;
 62    end
 63 endfunction
 64 
 65 
 66 function SymmetricTridiagonalQR(A,N)
 67    # QR for symmetric tridiagonal using givens rotations.  Will compute eigenpairs
 68    #  for matrices with random entries from 0 to 10.
 69    answer=1;
 70    A_initial=A;
 71    V0=eye(N);
 72    P0=A;
 73    m=1; # initialize m=number of iterations of while loop needed to reduce off-diagonals abs. value below chosen small number
 74    while (answer==1)
 75        for k=1:(N-1)
 76           Q0=eye(N);
 77           P=eye(N);
 78           [c,s]=givens_rotation(P0(k,k),P0(k+1,k));
 79           P(k,k)= c;
 80           P(k,k+1)= -s;
 81           P(k+1,k)= s;
 82           P(k+1,k+1)= c;
 83           P0=P*A;
 84           Q0=Q0*transpose(P);
 85           R0=P0;
 86           A=R0*Q0;
 87           V0=V0*Q0;
 88        endfor
 89        m=m+1;
 90        answer=CheckOffDiagaonal(A,N);
 91    endwhile
 92    for i=1:N
 93       lambda(i)=A(i,i);
 94    endfor
 95    for r=1:N
 96        # check eigenvalue equation equals approx. zero for each eigenpair
 97        eigenvalue_eq_check=(A_initial-(lambda(r)*eye(N)))*V0(:,r)
 98    endfor
 99 endfunction
100 
101 
102 function TestSymmetricTridiagonalQR(N)
103    A=GetA(N);
104    SymmetricTridiagonalQR(A,N);
105 endfunction
106 
107 
108 TestSymmetricTridiagonalQR(5);