octave:10> spdiags([[1;2;3;4;5],[6;7;8;9;10],[11;12;13;14;15]],[-1,0,1],5,5) ans = Compressed Column Sparse (rows = 5, cols = 5, nnz = 13 [52%]) (1, 1) -> 6 (2, 1) -> 1 (1, 2) -> 12 (2, 2) -> 7 (3, 2) -> 2 (2, 3) -> 13 (3, 3) -> 8 (4, 3) -> 3 (3, 4) -> 14 (4, 4) -> 9 (5, 4) -> 4 (4, 5) -> 15 (5, 5) -> 10 octave:11> full(spdiags([[1;2;3;4;5],[6;7;8;9;10],[11;12;13;14;15]],[-1,0,1],5,5)) ans = 6 12 0 0 0 1 7 13 0 0 0 2 8 14 0 0 0 3 9 15 0 0 0 4 10 octave:12> A=full(spdiags([[1;2;3;4;5],[6;7;8;9;10],[11;12;13;14;15]],[-1,0,1],5,5)) A = 6 12 0 0 0 1 7 13 0 0 0 2 8 14 0 0 0 3 9 15 0 0 0 4 10 octave:13> triu(A) ans = 6 12 0 0 0 0 7 13 0 0 0 0 8 14 0 0 0 0 9 15 0 0 0 0 10 octave:14> triu(A,1) ans = 0 12 0 0 0 0 0 13 0 0 0 0 0 14 0 0 0 0 0 15 0 0 0 0 0 octave:15> triu(A,1) ans = 0 12 0 0 0 0 0 13 0 0 0 0 0 14 0 0 0 0 0 15 0 0 0 0 0 octave:16> triu(A,0) ans = 6 12 0 0 0 0 7 13 0 0 0 0 8 14 0 0 0 0 9 15 0 0 0 0 10 octave:17> D=diag(diag(A)) D = Diagonal Matrix 6 0 0 0 0 0 7 0 0 0 0 0 8 0 0 0 0 0 9 0 0 0 0 0 10 octave:18> D=spdiags(diag(A),0,5,5) D = Compressed Column Sparse (rows = 5, cols = 5, nnz = 5 [20%]) (1, 1) -> 6 (2, 2) -> 7 (3, 3) -> 8 (4, 4) -> 9 (5, 5) -> 10 octave:19> D\A ans = 1.00000 2.00000 0.00000 0.00000 0.00000 0.14286 1.00000 1.85714 0.00000 0.00000 0.00000 0.25000 1.00000 1.75000 0.00000 0.00000 0.00000 0.33333 1.00000 1.66667 0.00000 0.00000 0.00000 0.40000 1.00000 octave:20> inv(D)*A % no ans = 1.00000 2.00000 0.00000 0.00000 0.00000 0.14286 1.00000 1.85714 0.00000 0.00000 0.00000 0.25000 1.00000 1.75000 0.00000 0.00000 0.00000 0.33333 1.00000 1.66667 0.00000 0.00000 0.00000 0.40000 1.00000 octave:21> A A = 6 12 0 0 0 1 7 13 0 0 0 2 8 14 0 0 0 3 9 15 0 0 0 4 10 octave:22> D D = Compressed Column Sparse (rows = 5, cols = 5, nnz = 5 [20%]) (1, 1) -> 6 (2, 2) -> 7 (3, 3) -> 8 (4, 4) -> 9 (5, 5) -> 10 octave:23> J=eye(5)-D\A J = 0.00000 -2.00000 -0.00000 -0.00000 -0.00000 -0.14286 0.00000 -1.85714 -0.00000 -0.00000 -0.00000 -0.25000 0.00000 -1.75000 -0.00000 -0.00000 -0.00000 -0.33333 0.00000 -1.66667 -0.00000 -0.00000 -0.00000 -0.40000 0.00000 octave:24> eig(J) ans = -1.2559e+00 -6.5012e-01 -1.0880e-16 1.2559e+00 6.5012e-01 octave:25> A = [2, 1, 0; 1, 2, 1; 0, 1, 2]; octave:26> b = A * ones (3, 1); octave:27> x = jacobi (A, b, 1e-8); octave:28> norm(x-ones(3,1)) ans = 1.8250e-08 octave:29> norm(x-ones(3,1))/norm(ones(3,1)) ans = 1.0537e-08 octave:30> [x,its] = jacobi (A, b, 1e-8); octave:31> its its = 53 octave:32> [x,its] = jacobi (A, b); octave:33> its its = 40 octave:34> diary off