QR method for eigenvalues

[nowy20180820]

Write your question here.

 // Example program
#include <iostream>
#include <iomanip>
#include <cmath>
#include <ctime>

const int maxiter = 1000;
const double eps = 1e-12;

void elimHess(int n,double **A)
{
    int i,j,l,k;
    double m,maxA,temp;
    for(i = 0;i < n - 2;i++)
    {
        k = i + 1;
        maxA = abs(A[k][i]);
        for(j = i + 1;j < n;j++)
           if(abs(A[j][i]) > maxA)
           {
               k = j;
               maxA = abs(A[k][i]);
           }
           if(maxA > 0)
           {
               for(j = 0;j < n;j++)
               {
                   temp = A[k][j];
                   A[k][j] = A[i + 1][j];
                   A[i + 1][j] = temp;
               }
               for(j = 0;j < n;j++)
               {
                   temp = A[j][k];
                   A[j][k] = A[j][i + 1];
                   A[j][i + 1] = temp;
               }
               for(j = i + 2;j < n;j++)
               {
                    m = (double)A[j][i]/A[i + 1][i];
                    for(l = 0; l < n;l++)
                       A[j][l] -= m * A[i + 1][l];
                    for(l = 0;l < n;l++)
                       A[l][i+1] += m * A[l][j];    
               }
           }
    }
}

void QRdecomp(int n,int m,double **A,double **Q)
{
    double r,c,s,temp;
    int i,j,k,min;
    for(i = 0;i < m;i++)
       for(j = 0;j < m;j++)
          Q[i][j] = (i == j ? 1 : 0);
    min = (m < n ? m : n);
    for(i = 0;i < min;i++)
        for(j = i + 1;j < m; j++)
           if(A[j][i] != 0)
           {
               r = hypot(A[i][i],A[j][i]);
               c = (double)A[i][i]/r;
               s = (double)A[j][i]/r;
               for(k = 0;k < n; k++)
               {
                   temp = A[i][k];
                   A[i][k] = c * A[i][k] + s * A[j][k];
                   A[j][k] = -s * temp + c * A[j][k];
               }
               for(k = 0 ;k < m;k++)
               {
                   temp = Q[k][i];
                   Q[k][i] = c * Q[k][i] + s * Q[k][j];
                   Q[k][j] = -s * temp + c * Q[k][j];
               }
           }
}

void multiplyMatrix(int n,int m ,int p,double **A,double **B,double **C)
{
    int i,j,k;
    for(i=0;i<n;i++)
        for(k = 0;k < p;k++)
        {
            C[i][k] = 0;
            for(j = 0;j < n;j++)
                C[i][k] += A[i][j] * B[j][k];
        }
}

void copyMatrix(int m,int n, double **A,double **B)
{
    for(int i=0; i< m;i++)
        for(int j=0;j < n;j++)
            B[i][j] = A[i][j];
}

void printMatrix(int m,int n,double**A)
{
    for(int i=0;i<n;i++)
    {
        for(int j=0;j<n;j++)
            if(abs(A[i][j]) > eps)
                std::cout<<std::setprecision(12)<<A[i][j]<<" ";
            else
                std::cout<<"0 ";
        std::cout<<"\n";
    }
    std::cout<<"\n";
}

int main()
{
    int m,n,k;
    double **A;
    double **Q;
    double **R;
    double p,q,d;
    std::cout<<"Podaj stopien macierzy \n";
    std::cin>>n;
    std::cout<<"Podaj gorny zakres przedzialu dla elementow macierzy \n";
    std::cin>>m;
    A = new double*[n];
    for(int i=0;i<n;i++)
       A[i] = new double[n];
    Q = new double*[n];
    for(int i=0;i<n;i++)
       Q[i] = new double[n];
    R = new double*[n];
    for(int i=0;i<n;i++)
       R[i] = new double[n];
    srand(time(NULL));
    for(int i=0;i<n;i++)
      for(int j=0;j<n;j++)
        A[i][j] = rand()%m;
    std::cout<<"Wylosowana macierz \n";
    printMatrix(n,n,A);
    elimHess(n,A);
    std::cout<<"Macierz po sprowadzeniu do postaci Hessenberga \n";
    printMatrix(n,n,A);
    for(int i=0;i<maxiter;i++)
    {
        QRdecomp(n,n,A,Q);
        copyMatrix(n,n,A,R);
        multiplyMatrix(n,n,n,R,Q,A);
    }
    std::cout<<"Macierz po iteracyjnym sprowadzeniu do postaci Schura \n";
    printMatrix(n,n,A);
    std::cout<<"Przyblizone wartosci wlasne macierzy \n";
    k = 0;
    while(k < n)
    {
        if(k + 1 < n && abs(A[k+1][k]) > eps)
        {
            p = 0.5 * (A[k][k] + A[k + 1][k + 1]);
            q = A[k][k] * A[k + 1][k + 1] - A[k][k + 1] * A[k+1][k];
            d = sqrt(abs(q - p * p));
            std::cout<<"x["<<k<<"]="<<std::setprecision(12)<<p<<"-"<<std::setprecision(12)<<d<<"*i \n";
            std::cout<<"x["<<k+1<<"]="<<std::setprecision(12)<<p<<"+"<<std::setprecision(12)<<d<<"*i \n";
            k += 2;
        }
        else
        {
            std::cout<<"x["<<k<<"]="<<std::setprecision(12)<<A[k][k]<<"\n";
            k++;
        }
    }
    for(int i=0;i<n;i++)
      delete[] A[i];
    delete[] A;
    for(int i=0;i<n;i++)
      delete[] Q[i];
    delete[] Q;
    for(int i=0;i<n;i++)
      delete[] R[i];
    delete[] R;
    return 0;
}

Here is my code for QR method for eigenvalues with reduction to Hessenberg form
It works but there some improvements can be made

I used Gaussian elimiation for Hessenberg reduction
There Householder reflections can be used but I don't know how to write code

1. Which shift should I choose
Doubly implicit shift is tempting because matrix has
real entries but eigenvalues can be complex
2. What are advantages of deflation
Can it be written in place (without extra space)
3. How can I accelerate convergence for repeated eigenvalues
Convergence is really slow in this case
4. What other stop criterion should I use

[jonnin]

I am really rusty but...
1) not sure right now
2) deflation is very important. its a double effect: you turn it into smaller problems (faster to solve) that can be solved independently (across multiple CPUs/threads, which is the name of the game since quad cores came out). So it will run faster than current time/#cpus
2.b) I think there are quick solves for 2x2 or 3x3 without all the iteration (?) so as your sub-divide, you can maybe break out earlier via a more direct approach?

3) Off the top of my head ... CHEAT :) If it 'looks like' it is converging on a previous value, stuff the known value in and move on. I don't know what criteria you need to decide though -- how close do they get before it starts bouncing around it? Can you know up front if there will be repeats? I seem to recall that you can maybe figure that out?

also, you can force zeros onto subdiagonals that are getting close. I can't tell if you already did that, though. This isn't just a blind assignment -- I think it requires a whole extra slew of code to do a transform and clear it out, maybe a givens transform (?). Sorry, its just been too long and there are too many approaches depending on what your data is.

Speaking of that, if this is for a specific subproblem set, you may have better choices. If its generic calculator, maybe not.

4) not sure right now.

[lastchance]

Even if your matrix itself is real I think you need to be working with complex<double>, as many (most?) of your eigenvalues and eigenvectors are going to be complex.