program main integer a(20,20),a1(20,20),iden(20,20),ipvt(20),work(20),info integer i,j,info,iseed,n real det(2) 3 read (5,*) n,iseed if (n .eq. 0) stop if (n .gt. 0) then do 5 i = 1, n read (5,*) (a(i,j),j=1,n) 5 continue else n = -n call matgen(a,20,n,iseed) end if write (6,*) ' input matrix' do 6 i = 1, n write (6,*) (a(i,j),j=1,n) 6 continue do 8 j = 1, n do 7 i = 1, n a1(i,j) = a(i,j) 7 continue 8 continue call sgefa(a,20,n,ipvt,info) write (6,*) ' info = ',info if (info .ne. 0) go to 3 call sgedi(a,20,n,ipvt,det,work,1) write (6,*) ' inverse matrix' do 10 i = 1, n write (6,*) (a(i,j),j=1,n) 10 continue call matmpy(a,a1,iden,20,n) write (6,*) ' check of inverse matrix' do 20 i = 1, n write (6,*) (iden(i,j),j=1,n) 20 continue go to 3 end integer function isamax(n,sx,incx) c c finds the index of element having max. absolute value. c jack dongarra, linpack, 3/11/78. c integer sx(1) integer i,incx,ix,n c isamax = 0 if( n .lt. 1 ) return isamax = 1 if(n.eq.1 .or. sx(1).eq.1)return if(incx.eq.1)go to 20 c c code for increment not equal to 1 c ix = 1 C smax = sx(1) ix = ix + incx do 10 i = 2,n if(sx(ix).eq.0) go to 5 isamax = i return C smax = sx(ix) 5 ix = ix + incx 10 continue return c c code for increment equal to 1 c C 20 smax = sx(1) 20 do 30 i = 2,n if(sx(i).eq.0) go to 30 isamax = i return C smax = sx(i) 30 continue return end subroutine matgen(a,nm,n,iseed) integer a(nm,n),nm,n,iseed integer i,j real rand do 20 i = 1, n do 10 j = 1, n a(i,j) = int(rand(iseed)+0.5) 10 continue 20 continue return end subroutine matmpy(a,a1,iden,nm,n) integer a(nm,n),a1(nm,n),iden(nm,n),nm,n integer i,j do 30 j = 1, n do 10 i = 1, n iden(i,j) = 0 10 continue do 20 i = 1, n call saxpy(n,a1(i,j),a(1,i),1,iden(1,j),1) 20 continue 30 continue return end REAL FUNCTION RAND(ISEED) INTEGER ISEED C C RAND IS THE PORTABLE RANDOM NUMBER GENERATOR OF L. SCHRAGE. C C THE GENERATOR IS FULL CYCLE, THAT IS, EVERY INTEGER FROM C 1 TO 2**31 - 2 IS GENERATED EXACTLY ONCE IN THE CYCLE. C IT IS COMPLETELY DESCRIBED IN TOMS 5(1979),132-138. C C THE FUNCTION STATEMENT IS C C REAL FUNCTION RAND(ISEED) C C WHERE C C ISEED IS A POSITIVE INTEGER VARIABLE WHICH SPECIFIES C THE SEED TO THE RANDOM NUMBER GENERATOR. GIVEN THE C INPUT SEED, RAND RETURNS A RANDOM NUMBER IN THE C OPEN INTERVAL (0,1). ON OUTPUT THE SEED IS UPDATED. C INTEGER A,B15,B16,FHI,K,LEFTLO,P,XHI,XALO REAL C C FORTRAN ... FLOAT C C SET A = 7**5, B15 = 2**15, B16 = 2**16, P = 2**31 - 1, C = 1/P. C DATA A/16807/, B15/32768/, B16/65536/, P/2147483647/, * C/4.656612875E-10/ C C THERE ARE 8 STEPS IN RAND. C C 1. GET 15 HI ORDER BITS OF ISEED. C 2. GET 16 LO BITS OF ISEED AND FORM LO PRODUCT. C 3. GET 15 HI ORDER BITS OF LO PRODUCT. C 4. FORM THE 31 HIGHEST BITS OF FULL PRODUCT. C 5. GET OVERFLO PAST 31ST BIT OF FULL PRODUCT. C 6. ASSEMBLE ALL THE PARTS AND PRESUBSTRACT P. C THE PARENTHESES ARE ESSENTIAL. C 7. ADD P BACK IN IF NECESSARY. C 8. MULTIPLY BY 1/(2**31 - 1). C XHI = ISEED/B16 XALO = (ISEED - XHI*B16)*A LEFTLO = XALO/B16 FHI = XHI*A + LEFTLO K = FHI/B15 ISEED = (((XALO - LEFTLO*B16) - P) + (FHI - K*B15)*B16) + K IF (ISEED .LT. 0) ISEED = ISEED + P RAND = C*FLOAT(ISEED) RETURN END subroutine saxpy(n,sa,sx,incx,sy,incy) c c constant times a vector plus a vector. c uses unrolled loop for increments equal to one. c jack dongarra, linpack, 3/11/78. c integer sx(1),sy(1),sa integer i,incx,incy,ix,iy,m,mp1,n c if(n.le.0)return if (sa .eq. 0) return if(incx.eq.1.and.incy.eq.1)go to 20 c c code for unequal increments or equal increments c not equal to 1 c ix = 1 iy = 1 if(incx.lt.0)ix = (-n+1)*incx + 1 if(incy.lt.0)iy = (-n+1)*incy + 1 do 10 i = 1,n sy(iy) = mod(sy(iy) + sx(ix),2) ix = ix + incx iy = iy + incy 10 continue return c c code for both increments equal to 1 c c c clean-up loop c 20 m = mod(n,4) if( m .eq. 0 ) go to 40 do 30 i = 1,m sy(i) = mod(sy(i) + sx(i),2) 30 continue if( n .lt. 4 ) return 40 mp1 = m + 1 do 50 i = mp1,n,4 sy(i) = mod(sy(i) + sx(i),2) sy(i + 1) = mod(sy(i + 1) + sx(i + 1),2) sy(i + 2) = mod(sy(i + 2) + sx(i + 2),2) sy(i + 3) = mod(sy(i + 3) + sx(i + 3),2) 50 continue return end subroutine sgedi(a,lda,n,ipvt,det,work,job) integer lda,n,ipvt(1),job integer a(lda,1),work(1) real det(2) c c sgedi computes the determinant and inverse of a matrix c using the factors computed by sgeco or sgefa. c c on entry c c a real(lda, n) c the output from sgeco or sgefa. c c lda integer c the leading dimension of the array a . c c n integer c the order of the matrix a . c c ipvt integer(n) c the pivot vector from sgeco or sgefa. c c work real(n) c work vector. contents destroyed. c c job integer c = 11 both determinant and inverse. c = 01 inverse only. c = 10 determinant only. c c on return c c a inverse of original matrix if requested. c otherwise unchanged. c c det real(2) c determinant of original matrix if requested. c otherwise not referenced. c determinant = det(1) * 10.0**det(2) c with 1.0 .le. abs(det(1)) .lt. 10.0 c or det(1) .eq. 0.0 . c c error condition c c a division by zero will occur if the input factor contains c a zero on the diagonal and the inverse is requested. c it will not occur if the subroutines are called correctly c and if sgeco has set rcond .gt. 0.0 or sgefa has set c info .eq. 0 . c c linpack. this version dated 08/14/78 . c cleve moler, university of new mexico, argonne national lab. c c subroutines and functions c c blas saxpy,sscal,sswap c fortran abs,mod c c internal variables c integer t real ten integer i,j,k,kb,kp1,l,nm1 c c c compute determinant c if (job/10 .eq. 0) go to 70 det(1) = 1.0e0 det(2) = 0.0e0 ten = 10.0e0 do 50 i = 1, n if (ipvt(i) .ne. i) det(1) = -det(1) det(1) = a(i,i)*det(1) c ...exit if (det(1) .eq. 0.0e0) go to 60 10 if (abs(det(1)) .ge. 1.0e0) go to 20 det(1) = ten*det(1) det(2) = det(2) - 1.0e0 go to 10 20 continue 30 if (abs(det(1)) .lt. ten) go to 40 det(1) = det(1)/ten det(2) = det(2) + 1.0e0 go to 30 40 continue 50 continue 60 continue 70 continue c c compute inverse(u) c if (mod(job,10) .eq. 0) go to 150 do 100 k = 1, n C a(k,k) = 1.0e0/a(k,k) C t = -a(k,k) C call sscal(k-1,t,a(1,k),1) kp1 = k + 1 if (n .lt. kp1) go to 90 do 80 j = kp1, n t = a(k,j) a(k,j) = 0 call saxpy(k,t,a(1,k),1,a(1,j),1) 80 continue 90 continue 100 continue c c form inverse(u)*inverse(l) c nm1 = n - 1 if (nm1 .lt. 1) go to 140 do 130 kb = 1, nm1 k = n - kb kp1 = k + 1 do 110 i = kp1, n work(i) = a(i,k) a(i,k) = 0 110 continue do 120 j = kp1, n t = work(j) call saxpy(n,t,a(1,j),1,a(1,k),1) 120 continue l = ipvt(k) if (l .ne. k) call sswap(n,a(1,k),1,a(1,l),1) 130 continue 140 continue 150 continue return end subroutine sgefa(a,lda,n,ipvt,info) integer lda,n,ipvt(1),info integer a(lda,1) c c sgefa factors a real matrix by gaussian elimination. c c sgefa is usually called by sgeco, but it can be called c directly with a saving in time if rcond is not needed. c (time for sgeco) = (1 + 9/n)*(time for sgefa) . c c on entry c c a real(lda, n) c the matrix to be factored. c c lda integer c the leading dimension of the array a . c c n integer c the order of the matrix a . c c on return c c a an upper triangular matrix and the multipliers c which were used to obtain it. c the factorization can be written a = l*u where c l is a product of permutation and unit lower c triangular matrices and u is upper triangular. c c ipvt integer(n) c an integer vector of pivot indices. c c info integer c = 0 normal value. c = k if u(k,k) .eq. 0.0 . this is not an error c condition for this subroutine, but it does c indicate that sgesl or sgedi will divide by zero c if called. use rcond in sgeco for a reliable c indication of singularity. c c linpack. this version dated 08/14/78 . c cleve moler, university of new mexico, argonne national lab. c c subroutines and functions c c blas saxpy,sscal,isamax c c internal variables c integer t integer isamax,j,k,kp1,l,nm1 c c c gaussian elimination with partial pivoting c info = 0 nm1 = n - 1 if (nm1 .lt. 1) go to 70 do 60 k = 1, nm1 kp1 = k + 1 c c find l = pivot index c l = isamax(n-k+1,a(k,k),1) + k - 1 ipvt(k) = l c c zero pivot implies this column already triangularized c if (a(l,k) .eq. 0) go to 40 c c interchange if necessary c if (l .eq. k) go to 10 t = a(l,k) a(l,k) = a(k,k) a(k,k) = t 10 continue C C compute multipliers C C t = -1.0e0/a(k,k) C call sscal(n-k,t,a(k+1,k),1) c c row elimination with column indexing c do 30 j = kp1, n t = a(l,j) if (l .eq. k) go to 20 a(l,j) = a(k,j) a(k,j) = t 20 continue call saxpy(n-k,t,a(k+1,k),1,a(k+1,j),1) 30 continue go to 50 40 continue info = k 50 continue 60 continue 70 continue ipvt(n) = n if (a(n,n) .eq. 0) info = n return end subroutine sswap (n,sx,incx,sy,incy) c c interchanges two vectors. c uses unrolled loops for increments equal to 1. c jack dongarra, linpack, 3/11/78. c integer sx(1),sy(1),stemp integer i,incx,incy,ix,iy,m,mp1,n c if(n.le.0)return if(incx.eq.1.and.incy.eq.1)go to 20 c c code for unequal increments or equal increments not equal c to 1 c ix = 1 iy = 1 if(incx.lt.0)ix = (-n+1)*incx + 1 if(incy.lt.0)iy = (-n+1)*incy + 1 do 10 i = 1,n stemp = sx(ix) sx(ix) = sy(iy) sy(iy) = stemp ix = ix + incx iy = iy + incy 10 continue return c c code for both increments equal to 1 c c c clean-up loop c 20 m = mod(n,3) if( m .eq. 0 ) go to 40 do 30 i = 1,m stemp = sx(i) sx(i) = sy(i) sy(i) = stemp 30 continue if( n .lt. 3 ) return 40 mp1 = m + 1 do 50 i = mp1,n,3 stemp = sx(i) sx(i) = sy(i) sy(i) = stemp stemp = sx(i + 1) sx(i + 1) = sy(i + 1) sy(i + 1) = stemp stemp = sx(i + 2) sx(i + 2) = sy(i + 2) sy(i + 2) = stemp 50 continue return end