LUdecomp_spec.F

c**********************************************************************
#include "author.inc"
c*    $Id: LUdecomp_spec.F,v 1.1 1995/11/14 02:19:42 turner Exp $
c*
c*    Subroutines which overwrite a square n by n matrix, a, stored in
c*    full format, by its LU decomposition.  KJI ordering of the loops
c*    is used.
c*
c*    The CVD$ compiler directives are for Alliant Fortran.
c*
c*    <PARAMETER LIST>
c*
c*     Input:
c*      idim - leading dimension of a
c*      n - number of unknowns
c*
c*     In/Out:
c*      a - matrix to be factored
c*
c*     Output:
c*      status - return status
c*        -4  ==>  breakdown
c*        -1  ==>  invalid argument(s)
c*         0  ==>  success
c*
c*    <SUBROUTINES REQUIRED>
c*
c*     JT_x_eq_sx
c*     JT_y_eq_y_plus_sx
c*
#include "copyright.inc"
c**********************************************************************
      subroutine JT_LUdecomp_Full (idim, n, a, status)
      implicit none
c
c ... Input:
      integer idim, n
c
c ... In/Out:
      real a(idim,n)
c
c ... Output:
      integer status
c
c ... Local:
      integer j, k
      real amult, one, zero
c
      parameter (zero=0.0d0, one=1.0d0)
c
c ... Check arguments.
      if (idim.le.0 .or. n.le.0 .or. n.gt.idim) then
       status = -1
       return
      endif
c
CVD$  NOCONCUR
CVD$  NOSYNC
      do k=1,n-1
       if (a(k,k) .eq. zero) then
        status = -4
        return
       endif
       amult = one / a(k,k)
       call JT_x_eq_sx (n-k, amult, a(k+1,k), status)
c
CVD$   NOCONCUR
       do j=k+1,n
        call JT_y_eq_y_plus_sx (n-k, -a(k,j), a(k+1,k), a(k+1,j),
     &       status)
       enddo
      enddo
c
      return
      end