Matrix Computation

2022 Fall

Benchmark problems

     The Koblenz Network Collection,        Stanford Large Network Dataset Collection

        E = load('butterfly.txt');

        n = max(max(E));

        m = size(E,2);

        if (m == 2 )

               mtx_L  = sparse(E(:,1), E(:,2), 1, n, n); 

        else

              mtx_L  = sparse(E(:,1), E(:,2), E(:,3), n, n);

        end

        D_L   = diag(mtx_L);

        if ( ~isempty(find(D_L ~= 0, 1)) )

             mtx_L = mtx_L - sparse(1:n,1:n,D_L);

       end        

       mtx_ L = mtx_L + mtx_L';

       D               = diag(sum(mtx_L, 2));

       mtx_L   = D - mtx_L;

      G                 = graph(E(:,1), E(:,2));

      [bins,binsizes]   = conncomp(G);

      idx_first_graph   = binsizes(bins) == max(binsizes); 

      mtx_L                       = mtx_L(idx_first_graph,idx_first_graph);

      SG = subgraph(G, idx_first_graph);

      plot(SG)

      %

       [V, D] = eig(full(mtx_L));

       trimesh(E,V(:,2),V(:,3),[],V(:,2),'LineWidth',3);

       hold on; pIDX = V(:,2)>0; nIDX = V(:,2)<0;

       plot(V(pIDX,2),V(pIDX,3),'ko','MarkerFaceColor','r');

       plot(V(nIDX,2),V(nIDX,3),'ko','MarkerFaceColor','b');

       view([0,0,1]); axis off;


Power and inverse power methods

            dim_n                                                  = size(mtx_L,1);

            sigma                                                     = 1.0e-3;

           mtx_shift_L                                      = mtx_L - sigma * speye(dim_n);

            mtx_factor.Perm_amd_vec     = amd(mtx_shift_L);

            W_reorder                                          = mtx_shift_L(mtx_factor.Perm_amd_vec, mtx_factor.Perm_amd_vec);

            mtx_factor.mtx_Perm_amd    = sparse(mtx_factor.Perm_amd_vec,1:dim_n,ones(dim_n,1));

            [mtx_factor.mtx_Low_L,mtx_factor.mtx_upper_U,mtx_factor.mtx_Perm_LU] = lu(W_reorder);

            sol = mtx_factor.mtx_Perm_amd * (mtx_factor.mtx_upper_U \ (mtx_factor.mtx_Low_L \ (mtx_factor.mtx_Perm_LU * rhs(mtx_factor.Perm_amd_vec,:))));

Krylov method for solving eigenvalue problem

     run_Lanczos_solver.m, Lanczos_basic.m, Lanczos_with_IR_Reorth.m, Lanczos_Steps_with_Reorth.m, Implicit_Restart_for_Lanczos.m


       ncv                   = 10;

        [EV, EW]       = eigs(@(x)Shift_A_inv_b(x,mtx_factor), dim_n, ncv, 'largestabs'); %,'SubspaceDimension',45);

        EW                    = sigma + 1./diag(EW);

        [~,idx_sort] = sort(abs(EW));

        EW                    = EW(idx_sort,1);

        EV                     = EV(:,idx_sort);

function [sol] = Shift_A_inv_b(rhs,mtx_factor)

sol = mtx_factor.mtx_Perm_amd * (mtx_factor.mtx_upper_U \ (mtx_factor.mtx_Low_L \ (mtx_factor.mtx_Perm_LU * rhs(mtx_factor.Perm_amd_vec,:))));

end

Conjugate gradient method

Preconditioning conjugate gradient method

Generalized Conjugate Gradient Methods for Non-symmetric Linear Systems

GMRES: Generalized Minimal Residual Algorithm for Solving Non-symmetric Linear Systems

Preconditioning

Polynomial Jacobi-Davidson Method

QR Algorithms