subroutine tgsld1(a1,a2,b,x,c1,c2,r1,r2,r3,r5,r6,r,m,l,lda) integer m,l,lda double precision a1(lda,l),a2(lda,1),b(m,l),x(m,l),c1(m,m,1), * c2(m,m,1),r1(m,m),r2(m,m),r3(m,m),r5(m,m),r6(m,m),r(m) c c tgsld1 solves the double precision linear system c a * x = b c with the tg - matrix a . c c on entry c c a1 double precision(m**2,l) c the first row of blocks of the tg - matrix a . c each block is represented by columns . c on return a1 is unaltered . c c a2 double precision(m**2,l - 1) c the first column of blocks of the tg - matrix a c beginning with the second block. each block is c represented by columns. on return a2 is unaltered . c c b double precision(m*l) c the right hand side vector . c on return b is unaltered . c c c1 double precision(m,m,l - 1) c a work array . c c c2 double precision(m,m,l - 1) c a work array . c c r1 double precision(m,m) c a work array . c c r2 double precision(m,m) c a work array . c c r3 double precision(m,m) c a work array . c c r5 double precision(m,m) c a work array . c c r6 double precision(m,m) c a work array . c c r double precision(m) c a work vector . c c m integer c the order of the blocks of the matrix a . c c l integer c the number of blocks in a row or column c of the matrix a . c c lda integer c the leading dimension of the array a . c c on return c c x double precision(m*l) c the solution vector. x may coincide with b . c c toeplitz package. this version dated 07/23/82 . c c subroutines and functions c c linpack ... daxpy,dgefa,dgesl c ... (for in-line daxpy, see directions in comments) c c internal variables c integer i,i1,i2,i3,ii,j,n,n1,n2 c c solve the system with the principal minor of order m . c i3 = 1 do 20 j = 1, m do 10 i = 1, m c1(i,j,1) = a1(i3,1) r1(i,j) = a1(i3,1) r3(i,j) = r1(i,j) i3 = i3 + 1 10 continue x(j,1) = b(j,1) 20 continue call dgefa(r3,m,m,r,ii) call dgesl(r3,m,m,r,x(1,1),0) if (l .eq. 1) go to 420 c c recurrent process for solving the system c with the tg - matrix for n = 2, l . c do 410 n = 2, l c c compute multiples of the first and last block columns of c the inverse of the principal minor of order m*n . c n1 = n - 1 n2 = n - 2 i3 = 1 do 40 j = 1, m do 30 i = 1, m r5(i,j) = a2(i3,n1) r6(i,j) = a1(i3,n) i3 = i3 + 1 30 continue 40 continue if (n .eq. 2) go to 100 do 60 j = 1, m do 50 i = 1, m c1(i,j,n1) = r2(i,j) 50 continue 60 continue do 90 i1 = 1, n2 i2 = n - i1 do 80 j = 1, m i3 = 1 do 70 i = 1, m c for in-line daxpy, activate next 5 lines and deactivate following 3 . c do 65 ii = 1, m c r5(ii,j) = r5(ii,j) + c1(i,j,i2)*a2(i3,i1) c r6(ii,j) = r6(ii,j) + c2(i,j,i1)*a1(i3,i1+1) c i3 = i3 + 1 c 65 continue call daxpy(m,c1(i,j,i2),a2(i3,i1),1,r5(1,j),1) call daxpy(m,c2(i,j,i1),a1(i3,i1+1),1,r6(1,j),1) i3 = i3 + m 70 continue 80 continue 90 continue 100 continue do 120 j = 1, m do 110 i = 1, m r2(i,j) = -r5(i,j) 110 continue call dgesl(r3,m,m,r,r2(1,j),0) 120 continue do 140 j = 1, m do 130 i = 1, m r3(i,j) = r6(i,j) r6(i,j) = -c1(i,j,1) 130 continue 140 continue do 160 j = 1, m do 150 i = 1, m c for in-line daxpy, activate next 3 lines and deactivate following 1 . c do 145 ii = 1, m c c1(ii,j,1) = c1(ii,j,1) + r2(i,j)*r3(ii,i) c 145 continue call daxpy(m,r2(i,j),r3(1,i),1,c1(1,j,1),1) 150 continue 160 continue call dgefa(r6,m,m,r,ii) do 180 j = 1, m call dgesl(r6,m,m,r,r3(1,j),0) do 170 i = 1, m c for in-line daxpy, activate next 3 lines and deactivate following 1 . c do 165 ii = 1, m c r1(ii,j) = r1(ii,j) + r3(i,j)*r5(ii,i) c 165 continue call daxpy(m,r3(i,j),r5(1,i),1,r1(1,j),1) 170 continue 180 continue if (n .eq. 2) go to 320 do 200 j = 1, m do 190 i = 1, m r6(i,j) = c2(i,j,1) 190 continue 200 continue do 310 i1 = 2, n1 if (i1 .eq. n1) go to 230 do 220 j = 1, m do 210 i = 1, m r5(i,j) = c2(i,j,i1) 210 continue 220 continue 230 continue do 260 j = 1, m do 240 i = 1, m c2(i,j,i1) = r6(i,j) 240 continue do 250 i = 1, m c for in-line daxpy, activate next 3 lines and deactivate following 1 . c do 245 ii = 1, m c c2(ii,j,i1) = c2(ii,j,i1) + r3(i,j)*c1(ii,i,i1) c 245 continue call daxpy(m,r3(i,j),c1(1,i,i1),1,c2(1,j,i1),1) 250 continue 260 continue do 280 j = 1, m do 270 i = 1, m c for in-line daxpy, activate next 3 lines and deactivate following 1 . c do 265 ii = 1, m c c1(ii,j,i1) = c1(ii,j,i1) + r2(i,j)*r6(ii,i) c 265 continue call daxpy(m,r2(i,j),r6(1,i),1,c1(1,j,i1),1) 270 continue 280 continue do 300 j = 1, m do 290 i = 1, m r6(i,j) = r5(i,j) 290 continue 300 continue 310 continue 320 continue do 340 j = 1, m do 330 i = 1, m c2(i,j,1) = r3(i,j) 330 continue 340 continue c c compute the solution of the system with the c principal minor of order m*n . c do 360 j = 1, m do 350 i = 1, m r3(i,j) = r1(i,j) 350 continue x(j,n) = b(j,n) 360 continue do 380 i1 = 1, n1 i2 = n - i1 i3 = 1 do 370 i = 1, m c for in-line daxpy, activate next 4 lines and deactivate following 2 . c do 365 ii = 1, m c x(ii,n) = x(ii,n) - x(i,i2)*a2(i3,i1) c i3 = i3 + 1 c 365 continue call daxpy(m,-x(i,i2),a2(i3,i1),1,x(1,n),1) i3 = i3 + m 370 continue 380 continue call dgefa(r3,m,m,r,ii) call dgesl(r3,m,m,r,x(1,n),0) do 400 i1 = 1, n1 do 390 i = 1, m c for in-line daxpy, activate next 3 lines and deactivate following 1 . c do 385 ii = 1, m c x(ii,i1) = x(ii,i1) + x(i,n)*c2(ii,i,i1) c 385 continue call daxpy(m,x(i,n),c2(1,i,i1),1,x(1,i1),1) 390 continue 400 continue 410 continue 420 continue return end