      subroutine nurbprocess(iobj)
       include 'ctparm.inc'
       include 'nrbparam.inc'
       include 'nurbsar.inc'
       include 'nrbmtrx.inc'
c
c initialize the flags that keep track of the 'current'
c nurbs matrices to 'invalid' so the first time 
c get_nurbs_matrices is called, a new set of matrices
c will be calculated
c
       nment = -1
       nmksegu = -1
       nmksegv = -1
c
      return
      end
c
       subroutine calc_nurbs_matrices(iobj,ksegu,ksegv)
c***************************************************************
c purpose: to compute the matrix product [mu][pw][mv]
c***************************************************************
       include 'ctparm.inc'
       include 'nrbparam.inc'
       include 'nurbsar.inc'
       include 'nrbmtrx.inc'
       dimension r(3)
       real*4 rwt(mxorder,mxorder),rwtcntrl(mxorder,mxorder,3),
     &        wtcntrl(mxcpts,mxcpts,3)
       real*8 tmp1(mxorder,mxorder),
     &      cbmtrxu(mxorder,mxorder),cbmtrxv(mxorder,mxorder),
     &      matrixu(mxkts,mxorder,mxcpts),
     &      matrixv(mxkts,mxorder,mxcpts)

c
      
      call blndmtrx(iorderuar(iobj),iulsuar(iobj),knotuar(1,iobj),
     & lstknotuar(iobj),matrixu)
      call blndmtrx(iordervar(iobj),iulsvar(iobj),knotvar(1,iobj),
     & lstknotvar(iobj),matrixv)
c
        call fndcbmfkseg(iorderuar(iobj),iulsuar(iobj),ksegu,
     &  matrixu,cbmtrxu)
        call fndcbmfkseg(iordervar(iobj),iulsvar(iobj),ksegv,
     &  matrixv,tmp1)
        call tpmtrxdd(mxorder,mxorder,tmp1,iordervar(iobj),
     &  iordervar(iobj),cbmtrxv)
       call curwtwtct(iobj,ksegu,ksegv,rwt,rwtcntrl)

c set up numerator
        do ixyz = 1,3
         call genmtrxmuldsd(mxorder,mxorder,mxorder,cbmtrxu,
     &   rwtcntrl(1,1,ixyz),iorderuar(iobj),iorderuar(iobj),
     &   iordervar(iobj),tmp1)
         call genmtrxmuldds(mxorder,mxorder,mxorder,tmp1,cbmtrxv,
     &   iorderuar(iobj),iordervar(iobj),iordervar(iobj),
     &   bwcpuvm(1,1,ixyz))
        enddo

c set up denominator
         call genmtrxmuldsd(mxorder,mxorder,mxorder,cbmtrxu,rwt,
     &   iorderuar(iobj),iorderuar(iobj),iordervar(iobj),tmp1)
         call genmtrxmuldds(mxorder,mxorder,mxorder,tmp1,cbmtrxv,
     &   iorderuar(iobj),iordervar(iobj),iordervar(iobj),
     &   bwuvm(1,1))
c
      return
      end
c
      subroutine get_nurbs_matrices(iobj,ksegu,ksegv)
      include 'ctparm.inc'
      include 'nrbparam.inc'
      include 'nurbsar.inc'
      include 'nrbmtrx.inc'
c
      integer iobj,ksegu,ksegv
c      
      if( (iobj.ne.nment).or.(ksegu.ne.nmksegu).or.
     &(ksegv.ne.nmksegv) ) then
        call calc_nurbs_matrices(iobj,ksegu,ksegv)
      endif

      nment = iobj
      nmksegu = ksegu
      nmksegv = ksegv
c
      return
      end
c
       subroutine getvertex(iobj,u,v,r)
c***************************************************************
c purpose: to compute the catresian coordinates of a point (u,v) 
c on a nurb surface. 
c***************************************************************
       include 'ctparm.inc'
       include 'nrbparam.inc'
       include 'nurbsar.inc'
       include 'nrbmtrx.inc'
       dimension r(3)
       real*8 tmp1(mxorder,mxorder),
     &      num(mxorder,mxorder),denom(mxorder,mxorder)


c get appropriate matrices
       call getksegfuv(mxkts,knotuar(1,iobj),lstknotuar(iobj),u,ksegu)
       call getksegfuv(mxkts,knotvar(1,iobj),lstknotvar(iobj),v,ksegv)
       call get_nurbs_matrices(iobj,ksegu,ksegv)
c set up denominator
       call mulu(mxorder,u,bwuvm(1,1),iorderuar(iobj),
     & iordervar(iobj),tmp1)
       call mulv(mxorder,tmp1,1,iordervar(iobj),v,denom)
c set up numerator
c perform division and compute x, y and z
c (note : in the nurb expression 0/0=0)
c
       do ixyz = 1,3
         r(ixyz) = 0
         call mulu(mxorder,u,bwcpuvm(1,1,ixyz),
     &   iorderuar(iobj),iordervar(iobj),tmp1)
         call mulv(mxorder,tmp1,1,iordervar(iobj),v,num)
         if(denom(1,1).ne.0) r(ixyz) = num(1,1)/denom(1,1)
      enddo
c
c      write(6,*) 'iobj,u,v,r',iobj,u,v,r(1),r(2),r(3)

       return
       end

       subroutine gettans(iobj,u,v,r,ru,rv)
c****************************************************************
c purpose: to compute the catresian coordinates of a point (u,v)
c	   and its two tangent vectors on a nurb surface.
c****************************************************************

       include 'ctparm.inc'
       include 'nrbparam.inc'
       include 'nurbsar.inc'
       include 'nrbmtrx.inc'
       dimension r(3),ru(3),rv(3)
       real*8 tmp1(mxorder,mxorder),
     &      num(mxorder,mxorder),denom(mxorder,mxorder),
     &      numdu(mxorder,mxorder),dendu(mxorder,mxorder),
     &      numdv(mxorder,mxorder),dendv(mxorder,mxorder),
     &      numduv(mxorder,mxorder),denduv(mxorder,mxorder)

c

c get appropriate matrices
       call getksegfuv(mxkts,knotuar(1,iobj),lstknotuar(iobj),u,ksegu)
       call getksegfuv(mxkts,knotvar(1,iobj),lstknotvar(iobj),v,ksegv)
       call get_nurbs_matrices(iobj,ksegu,ksegv)
c set up denominator
       call mulu(mxorder,u,bwuvm(1,1),iorderuar(iobj),
     & iordervar(iobj),tmp1)
       call mulv(mxorder,tmp1,1,iordervar(iobj),v,denom)
       call mulvdv(tmp1,1,iordervar(iobj),v,dendv)
       call muludu(u,bwuvm(1,1),iorderuar(iobj),
     & iordervar(iobj),tmp1)
       call mulv(mxorder,tmp1,1,iordervar(iobj),v,dendu)

c set up numerator
c perform division and compute x, y and z
c the derivatives dr/du & dr/dv are computed by the derivative of
c the quotient  d(n/d)/du(v)
c (note : in the nurb expression 0/0=0)
c
       do k = 1,3
         r(k) = 0.	!x,y,z
         ru(k) = 0. 
         rv(k) = 0.
         call mulu(mxorder,u,bwcpuvm(1,1,k),
     &   iorderuar(iobj),iordervar(iobj),tmp1)
         call mulv(mxorder,tmp1,1,iordervar(iobj),v,num)
         call mulvdv(tmp1,1,iordervar(iobj),v,numdv)
         call muludu(u,bwcpuvm(1,1,k),iorderuar(iobj),
     &   iordervar(iobj),tmp1)
         call mulv(mxorder,tmp1,1,iordervar(iobj),v,numdu)
         if(denom(1,1).ne.0) then
               xden = denom(1,1)
               xnum = num(1,1)
               xdenu = dendu(1,1)
               xnumu = numdu(1,1)
               xdenv = dendv(1,1)
               xnumv = numdv(1,1)
               r(k) = xnum/xden
               ru(k) = ( xnumu - r(k) * xdenu ) / xden
               rv(k) = ( xnumv - r(k) * xdenv ) / xden
       endif
      enddo
c
       return
       end
c
       subroutine gethotans(iobj,u,v,r,ru,rv,ruu,rvv,ruv)
c****************************************************************
c purpose: to compute the catresian coordinates of a point (u,v),
c	   its two tangent vectors, and three 2nd order 
c          derivatives on a nurb surface.
c****************************************************************

       include 'ctparm.inc'
       include 'nrbparam.inc'
       include 'nurbsar.inc'
       include 'nrbmtrx.inc'
       dimension r(3),ru(3),rv(3),ruu(3),rvv(3),ruv(3)
       real*8 tmp1(mxorder,mxorder),
     &      num(mxorder,mxorder),denom(mxorder,mxorder),
     &      numdu(mxorder,mxorder),dendu(mxorder,mxorder),
     &      numdv(mxorder,mxorder),dendv(mxorder,mxorder),
     &      numduv(mxorder,mxorder),denduv(mxorder,mxorder),
     &      numduu(mxorder,mxorder),denduu(mxorder,mxorder),
     &      numdvv(mxorder,mxorder),dendvv(mxorder,mxorder)

c

c get appropriate matrices
       call getksegfuv(mxkts,knotuar(1,iobj),lstknotuar(iobj),u,ksegu)
       call getksegfuv(mxkts,knotvar(1,iobj),lstknotvar(iobj),v,ksegv)
       call get_nurbs_matrices(iobj,ksegu,ksegv)
c set up denominator
       call mulu(mxorder,u,bwuvm(1,1),iorderuar(iobj),
     & iordervar(iobj),tmp1)
       call mulv(mxorder,tmp1,1,iordervar(iobj),v,denom)
       call mulvdv(tmp1,1,iordervar(iobj),v,dendv)
       call mulvdvv(tmp1,1,iordervar(iobj),v,dendvv)
       call muludu(u,bwuvm(1,1),iorderuar(iobj),
     & iordervar(iobj),tmp1)
       call mulv(mxorder,tmp1,1,iordervar(iobj),v,dendu)
       call mulvdv(tmp1,1,iordervar(iobj),v,denduv)
       call muluduu(u,bwuvm(1,1),iorderuar(iobj),
     & iordervar(iobj),tmp1)
       call mulv(mxorder,tmp1,1,iordervar(iobj),v,denduu)

c set up numerator
c perform division and compute x, y and z
c the derivatives dr/du & dr/dv are computed by the derivative of
c the quotient  d(n/d)/du(v)
c (note : in the nurb expression 0/0=0)
c
       do k = 1,3
         r(k) = 0.	!x,y,z
         ru(k) = 0. 
         rv(k) = 0.
         ruu(k) = 0.
         rvv(k) = 0.
         ruv(k) = 0.
         call mulu(mxorder,u,bwcpuvm(1,1,k),
     &   iorderuar(iobj),iordervar(iobj),tmp1)
         call mulv(mxorder,tmp1,1,iordervar(iobj),v,num)
         call mulvdv(tmp1,1,iordervar(iobj),v,numdv)
         call mulvdvv(tmp1,1,iordervar(iobj),v,numdvv)
         call muludu(u,bwcpuvm(1,1,k),iorderuar(iobj),
     &   iordervar(iobj),tmp1)
         call mulv(mxorder,tmp1,1,iordervar(iobj),v,numdu)
         call mulvdv(tmp1,1,iordervar(iobj),v,numduv)
         call muluduu(u,bwcpuvm(1,1,k),iorderuar(iobj),
     &   iordervar(iobj),tmp1)
         call mulv(mxorder,tmp1,1,iordervar(iobj),v,numduu)

         if(denom(1,1).ne.0) then
               xden = denom(1,1)
               xnum = num(1,1)
               xdenu = dendu(1,1)
               xnumu = numdu(1,1)
               xdenv = dendv(1,1)
               xnumv = numdv(1,1)
               xdenuu = denduu(1,1)
               xnumuu = numduu(1,1)
               xdenvv = dendvv(1,1)
               xnumvv = numdvv(1,1)
               xdenuv = denduv(1,1)
               xnumuv = numduv(1,1)
               r(k) = xnum/xden
               ru(k) = ( xnumu - r(k) * xdenu ) / xden
               rv(k) = ( xnumv - r(k) * xdenv ) / xden
               ruu(k) =( xnumuu - 2.*ru(k)*xdenu - r(k)*xdenuu ) / xden
               rvv(k) =( xnumvv - 2.*rv(k)*xdenv - r(k)*xdenvv ) / xden
               ruv(k) =( xnumuv - rv(k)*xdenu - ru(k)*xdenv -  
     &		                  r(k)*xdenuv ) / xden
       endif
      enddo
c
       return
       end
c
      subroutine mulu(mx,u,matrix2,irow2,icolumn2,prod)
      real*8 matrix2(mx,mx),prod(mx,mx)
c
      zep = 1e-06
c multiply a given matrix by the u (a single row) matrix
      do icolumn = 1,icolumn2
        prod(1,icolumn) = 0.0
      enddo
      do icolumn = 1,icolumn2
         do irow = 1,irow2
            ipwr = irow2-irow
            if(u .eq. 0.0 .and. ipwr .eq. 0) then 
              prod(1,icolumn) = prod(1,icolumn)+ 
     &                               matrix2(irow,icolumn)
            else
              prod(1,icolumn) = prod(1,icolumn)+ 
     &                     matrix2(irow,icolumn)*(u**ipwr)
            endif
         enddo
         if(abs(prod(1,icolumn)).le.zep) prod(1,icolumn) = 0.0
      enddo
      return
      end
c
      subroutine mulv(mx,matrix1,irow1,icolumn1,v,prod)
      real*8 matrix1(mx,mx),prod(mx,mx)
c
      zep = 1e-06
c multiply a given matrix by the v (a single column) matrix
      do irow = 1,irow1
         prod(irow,1) = 0.0
      enddo
      do irow = 1,irow1
         do icolumn = 1,icolumn1
            ipwr = icolumn1-icolumn
            if(v .eq. 0.0 .and. ipwr .eq. 0) then 
               prod(irow,1) = prod(irow,1) + 
     &                    matrix1(irow,icolumn)
            else
               prod(irow,1) = prod(irow,1) + 
     &          matrix1(irow,icolumn)*(v**ipwr)
            endif
         enddo
         if(abs(prod(irow,1)).le.zep) prod(irow,1) = 0.0
      enddo
      return
      end

      subroutine muludu(u,matrix2,irow2,icolumn2,prod)

      include 'nrbparam.inc'
      real*8 matrix2(mxorder,mxorder),prod(mxorder,mxorder)
c
      zep = 1e-06
c multiply a given matrix by the derivative of the 
c u (a single row) matrix
      do icolumn = 1,icolumn2
        prod(1,icolumn) = 0.0
      enddo

      do icolumn = 1,icolumn2
         do irow = 1,irow2-1		!(last column of [u]'=0)
            ipwr = irow2-irow
            ipwr1 = ipwr-1
               if(u .eq. 0.0 .and. ipwr1 .eq. 0) then 
                  prod(1,icolumn) = prod(1,icolumn)+ 
     &                             matrix2(irow,icolumn)
               else
                  prod(1,icolumn) = prod(1,icolumn)+ 
     &                   matrix2(irow,icolumn)*(ipwr*u**ipwr1)
               endif
         enddo
         if(abs(prod(1,icolumn)).le.zep) prod(1,icolumn) = 0.0
      enddo
      return
      end
c
      subroutine mulvdv(matrix1,irow1,icolumn1,v,prod)

      include 'nrbparam.inc'
      real*8 matrix1(mxorder,mxorder),prod(mxorder,mxorder)
c
      zep = 1e-06
c multiply a given matrix by the derivative of the 
c v (a single column) matrix
      do irow = 1,irow1
        prod(irow,1) = 0.0
      enddo
      do irow = 1,irow1
         do icolumn = 1,icolumn1-1	!(last row of v' = 0)
            ipwr = icolumn1-icolumn
            ipwr1 = ipwr-1
               if(v .eq. 0.0 .and. ipwr1 .eq. 0) then 
                  prod(irow,1) = prod(irow,1) + 
     &                       matrix1(irow,icolumn)
               else
                  prod(irow,1) = prod(irow,1) + 
     &             matrix1(irow,icolumn)*(ipwr*v**ipwr1)
               endif
         enddo
         if(abs(prod(irow,1)).le.zep) prod(irow,1) = 0.0
      enddo
      return
      end


      subroutine muluduu(u,matrix2,irow2,icolumn2,prod)

      include 'nrbparam.inc'
      real*8 matrix2(mxorder,mxorder),prod(mxorder,mxorder)
c
      zep = 1e-06
c multiply a given matrix by the 2nd derivative of the 
c u (a single row) matrix
      do icolumn = 1,icolumn2
        prod(1,icolumn) = 0.0
      enddo

      do icolumn = 1,icolumn2
         do irow = 1,irow2-2		!(last 2 columns of [u]''=0)
            ipwr = irow2-irow
            ipwr1 = ipwr-1
            ipwr2 = ipwr-2
               if(u .eq. 0.0 .and. ipwr2 .eq. 0) then 
                  prod(1,icolumn) = prod(1,icolumn)+ 
     &                             matrix2(irow,icolumn) * 2.
               else
                  prod(1,icolumn) = prod(1,icolumn)+ 
     &                       matrix2(irow,icolumn) * 
     &                      (ipwr*ipwr1*u**ipwr2)
               endif
         enddo
         if(abs(prod(1,icolumn)).le.zep) prod(1,icolumn) = 0.0
      enddo
      return
      end
c
      subroutine mulvdvv(matrix1,irow1,icolumn1,v,prod)

      include 'nrbparam.inc'
      real*8 matrix1(mxorder,mxorder),prod(mxorder,mxorder)
c
      zep = 1e-06
c multiply a given matrix by the 2nd derivative of the 
c v (a single column) matrix
      do irow = 1,irow1
        prod(irow,1) = 0.0
      enddo
      do irow = 1,irow1
         do icolumn = 1,icolumn1-2	!(last 2 rows of [v]'' = 0)
            ipwr = icolumn1-icolumn
            ipwr1 = ipwr-1
            ipwr2 = ipwr-2
               if(v .eq. 0.0 .and. ipwr2 .eq. 0) then 
                  prod(irow,1) = prod(irow,1) + 
     &                       matrix1(irow,icolumn) * 2.
               else
                  prod(irow,1) = prod(irow,1) + 
     &                      matrix1(irow,icolumn)*
     &                     (ipwr*ipwr1*v**ipwr2)
               endif
         enddo
         if(abs(prod(irow,1)).le.zep) prod(irow,1) = 0.0
      enddo
      return
      end

      subroutine blndmtrx(iorder,iuls,knot,lstknot,matrix)
c*********************************************************
      include 'nrbparam.inc'
      real*8 coeff(mxkts,mxorder,mxorder),
     &     matrix(mxkts,mxorder,mxcpts),
     &     a(mxorder),b(mxorder)
      real*4 knot(mxkts)
c
      zep = 1e-06
      do imatnum = 1,lstknot - 1
         irow = iorder + 1
c
         do ii = 1,lstknot-1
           do ipwr1 = 1,iorder
             if (ii .eq. imatnum .and. knot(ii) .ne. knot(ii+1)) then
                coeff(ii,1,ipwr1) = 1
             else 
                coeff(ii,1,1) = 0.0
             endif
           enddo
         enddo
c
         do k = 2,iorder
            do ipwr1 = 1,k
               do i = 0,lstknot-1
               ii = i + 1
               a(k) = knot(ii+k-1) - knot(ii)
               b(k) = knot(ii+k) - knot(ii+1)
               if(abs(a(k)).le.zep) a(k) = 0.0
               if(abs(b(k)).le.zep) b(k) = 0.0
               if(ipwr1 .eq. 1) then
                 coeff(ii,k,ipwr1) = 0.0
                 if (b(k) .ne. 0.0) 
     &                coeff(ii,k,ipwr1) = ( knot(ii+k) * 
     &                        coeff(ii+1,k-1,1) )/b(k) 
                 if (a(k) .ne. 0.0) coeff(ii,k,ipwr1) = 
     &              coeff(ii,k,ipwr1)-knot(ii)*coeff(ii,k-1,1)/a(k)
               elseif(ipwr1 .eq. k) then
                 coeff(ii,k,ipwr1) = 0.0
                 if(a(k) .ne. 0.0)
     &              coeff(ii,k,ipwr1) = coeff(ii,k-1,k-1)/a(k)
                 if(b(k) .ne. 0.0) coeff(ii,k,ipwr1) =
     &              coeff(ii,k,ipwr1) - coeff(ii+1,k-1,k-1)/b(k)
               else
                  coeff(ii,k,ipwr1) = 0.0
                  if(a(k) .ne. 0.0)
     &              coeff(ii,k,ipwr1) = coeff(ii,k-1,ipwr1-1)/a(k) -
     &                           knot(ii)*coeff(ii,k-1,ipwr1)/a(k)
                  if(b(k) .ne. 0.0) coeff(ii,k,ipwr1) =
     &              coeff(ii,k,ipwr1) - coeff(ii+1,k-1,ipwr1-1)/b(k) +
     &                         knot(ii+k)*coeff(ii+1,k-1,ipwr1)/b(k) 
               endif 
               enddo
               if(k .eq. iorder) then
                  irow = irow - 1
                  do icolumn = 1,iuls
                    matrix(imatnum,irow,icolumn)=coeff(icolumn,k,ipwr1)
                  enddo
               endif
            enddo
         enddo
      enddo
      return
      end


      subroutine fndcbmtrx(u,iorder,iuls,lstknot,knot,blndmtrx,
     & knotseg,cbmtrx)
      include 'nrbparam.inc'
      real*8 blndmtrx(mxkts,mxorder,mxcpts),
     &     cbmtrx(mxorder,mxorder),knot(mxkts)
c
c   determine which knot segment u belongs to
      call getksegfuv(mxkts,knot,lstknot,u,knotseg)
c
c determine current blending function matrix
      call fndcbmfkseg(iorder,iuls,knotseg,blndmtrx,cbmtrx)
c
      return
      end
c
      subroutine getksegfuv(mxkts,knot,lstknot,uv,knotseg)
      real knot(mxkts)
c
         do kseg = 1,lstknot-1
            if ( knot(kseg+1) .ne. knot(kseg) ) then
               if ( uv .ge. knot(kseg) .and.
     &           uv .le. knot(kseg+1) ) knotseg = kseg
            endif
         enddo
c
      return
      end
c
       subroutine curwtwtct(iobj,ksegu,ksegv,rwt,rwtcntrl)
       include 'ctparm.inc'
       include 'nrbparam.inc'
       include 'nurbsar.inc'
       real rwt(mxorder,mxorder),rwtcntrl(mxorder,mxorder,3),
     & wtcntrl(mxcpts,mxcpts,3)
c
c determine reduced wt and wtcntrl matrices 
c
c generate matrix of weighted control points
      do irow = 1,iulsuar(iobj)
         do icolumn = 1,iulsvar(iobj)
            do ixyz = 1,3
            wtcntrl(irow,icolumn,ixyz) = 
     &      wtar(irow,icolumn,iobj) * cntrlptar(irow,icolumn,ixyz,iobj)
            enddo
         enddo
      enddo
c
           do i = 1,iorderuar(iobj)
              do j =1,iordervar(iobj)
                 ii = i+ksegu-iorderuar(iobj)
                 jj = j+ksegv-iordervar(iobj)
                 rwt(i,j)=wtar(ii,jj,iobj)
                 do ixyz = 1,3
                    rwtcntrl(i,j,ixyz)=wtcntrl(ii,jj,ixyz)
                 enddo
              enddo
           enddo

       return
       end

      subroutine fndcbmfkseg(iorder,iuls,kseg,blndmtrx,cbmtrx)
      include 'nrbparam.inc'
      real*8 blndmtrx(mxkts,mxorder,mxcpts),
     &     cbmtrx(mxorder,mxorder),knot(mxkts)
c
c determine current blending function matrix
        do irow = 1,iorder
           do icolumn = 1,iorder
              i = irow
              j = kseg-iorder+icolumn
              cbmtrx(irow,icolumn) = blndmtrx(kseg,i,j)
           enddo
        enddo
c
      return
      end
c
      subroutine genmtrxmul(mxr1,mxc1,mxc2,matrix1,matrix2,
     &                      irow1,icolumn1,icolumn2,prod)
      real matrix1(mxr1,mxc1),matrix2(mxc1,mxc2),
     &     prod(mxr1,mxc2)
c
       do irow = 1,irow1
           do icolumn = 1,icolumn2
             prod(irow,icolumn) = 0.0
           enddo
       enddo
c
       do irow =1,irow1
          do icolumn = 1,icolumn2
             do icolrow =1,icolumn1
                prod(irow,icolumn) = prod(irow,icolumn)+
     &                            matrix1(irow,icolrow)*
     &                         matrix2(icolrow,icolumn)
            enddo
         enddo
      enddo
      return
      end
c
      subroutine genmtrxmuldsd(mxr1,mxc1,mxc2,matrix1,matrix2,
     &                      irow1,icolumn1,icolumn2,prod)
      real*8 matrix1(mxr1,mxc1),prod(mxr1,mxc2)
      real*4 matrix2(mxc1,mxc2)
c
       do irow = 1,irow1
           do icolumn = 1,icolumn2
             prod(irow,icolumn) = 0.0
           enddo
       enddo
c
       do irow =1,irow1
          do icolumn = 1,icolumn2
             do icolrow =1,icolumn1
                prod(irow,icolumn) = prod(irow,icolumn)+
     &                            matrix1(irow,icolrow)*
     &                         matrix2(icolrow,icolumn)
            enddo
         enddo
      enddo
      return
      end
c
      subroutine genmtrxmuldds(mxr1,mxc1,mxc2,matrix1,matrix2,
     &                      irow1,icolumn1,icolumn2,prod)
      real*8 matrix1(mxr1,mxc1),matrix2(mxc1,mxc2)
      real*8 prod(mxr1,mxc2)
c
       do irow = 1,irow1
           do icolumn = 1,icolumn2
             prod(irow,icolumn) = 0.0
           enddo
       enddo
c
       do irow =1,irow1
          do icolumn = 1,icolumn2
             do icolrow =1,icolumn1
                prod(irow,icolumn) = prod(irow,icolumn)+
     &                            matrix1(irow,icolrow)*
     &                         matrix2(icolrow,icolumn)
            enddo
         enddo
      enddo
      return
      end
c
      subroutine tpmtrx(mxr,mxc,amtrx,lstrow,lstcol,transamtrx)
      real amtrx(mxr,mxc),transamtrx(mxr,mxc)
      do irow = 1,lstcol
         do icol = 1,lstrow
            transamtrx(irow,icol) = amtrx(icol,irow)
         enddo
      enddo
      return
      end
c
      subroutine tpmtrxdd(mxr,mxc,amtrx,lstrow,lstcol,transamtrx)
      real*8 amtrx(mxr,mxc),transamtrx(mxr,mxc)
      do irow = 1,lstcol
         do icol = 1,lstrow
            transamtrx(irow,icol) = amtrx(icol,irow)
         enddo
      enddo
      return
      end

      subroutine nurbsar(ibox)
      include 'ctparm.inc'
      include 'nrbparam.inc'
      include 'nurbsar.inc'
      include 'nrbcnst.inc'
      include 'nrbar.inc'
c
      ubegar(ibox) = ubeg
      uendar(ibox) = uend
      vbegar(ibox) = vbeg
      vendar(ibox) = vend
      iorderuar(ibox) = iorderu
      iordervar(ibox) = iorderv
      iulsuar(ibox) = iulsu
      iulsvar(ibox) = iulsv
      idegreeuar(ibox) = idegreeu
      idegreevar(ibox) = idegreev
      icluar(ibox) = iclu
      iclvar(ibox) = iclv
      iperuar(ibox) = iperu
      ipervar(ibox) = iperv
      lstknotuar(ibox) = lstknotu
      lstknotvar(ibox) = lstknotv
c
      do i = 1,lstknotu
         knotuar(i,ibox) = knotu(i)
      enddo
      do i = 1,lstknotv
         knotvar(i,ibox) = knotv(i)
      enddo
c
      do i = 1,iulsu
         do j = 1,iulsv
            wtar(i,j,ibox) = wt(i,j)
            do ixyz = 1,3
               cntrlptar(i,j,ixyz,ibox) = cntrlpt(i,j,ixyz)
c               write(6,*) 'ibox,cpt',ibox,cntrlptar(i,j,ixyz,ibox)
            enddo
c            write(6,*)
         enddo
      enddo
c
      return
      end
      subroutine nurbars(ibox)
      include 'nrbparam.inc'
      include 'ctparm.inc'
      include 'nurbsar.inc'
      include 'nrbcnst.inc'
      include 'nrbar.inc'
c
      ubeg = ubegar(ibox) 
      uend = uendar(ibox) 
      vbeg = vbegar(ibox) 
      vend = vendar(ibox) 
      iorderu = iorderuar(ibox) 
      iorderv = iordervar(ibox) 
      iulsu = iulsuar(ibox) 
      iulsv = iulsvar(ibox) 
      idegreeu = idegreeuar(ibox) 
      idegreev = idegreevar(ibox) 
      iclu = icluar(ibox) 
      iclv = iclvar(ibox) 
      iperu = iperuar(ibox) 
      iperv = ipervar(ibox) 
      lstknotu = lstknotuar(ibox) 
      lstknotv = lstknotvar(ibox) 
c
      do i = 1,lstknotu
         knotu(i) = knotuar(i,ibox) 
      enddo
      do i = 1,lstknotv
         knotv(i) = knotvar(i,ibox) 
      enddo
c
      do i = 1,iulsu
         do j = 1,iulsv
            wt(i,j) = wtar(i,j,ibox)
            do ixyz = 1,3
               cntrlpt(i,j,ixyz) = cntrlptar(i,j,ixyz,ibox)
c               write(6,*) 'ibox,cptsc',ibox,cntrlpt(i,j,ixyz)
            enddo
c            write(6,*)
         enddo
      enddo
c
      return
      end
