program puz16
!    use LAPACK, only:sgesv
    parameter upr=16,n=3,nrhs=4,lda=3,ldb=4
    integer info,ipiv(99)
    integer ok(16)/16*0/
    integer nrow(4)/4*0/
    integer ntry/0/
    real a(3,3),b(3)
    CALL SOLVE(A,3,3,B)
    do i1=1,upr
        ok(1)=i1
        do i2=1,upr
            ok(2)=i2
            if(i2.ne.i1)then
           
            do i3=1,upr
                ok(3)=i3
                if(i3.ne.i2 .and. i3.ne.i1)then
                i4=36-i1-i2-i3
                if(chek(3,i4) .and. i4>0 .and. i4 <17)then
                    ok(4)=i4
                    print 104,ok(1:4)
104                 format(4i4)
                    nrow(1)=nrow(1)+1
                    call row2()
                endif
                endif
                
            enddo !i3
            endif
          
       enddo ! i2
       
    enddo !i1
    print *,"nrow=",nrow
    read(*,*)
    contains
    logical function chek(mm,nn)
    integer ii,mm,nn
    chek=.true.
    do ii=1,mm
        if(ok(ii)==nn)go to 90
    enddo
    return
90  chek=.false.
    end function chek
    subroutine row2()
    integer i1,i2,i3,i4
    do i1=1,upr
        ok(5)=i1
        if(chek(4,i1))then
        
        do i2=1,upr
            ok(6)=i2
            if(chek(5,i2))then
           
            do i3=1,upr                
                ok(7)=i3
                if(chek(6,i3))then
        ! complete row 3
                    i4=36-i1-i2-i3                   
                    if(chek(7,I4).and. &
                      i4>0 .and. i4 <17)then
                    ok(8)=i4
                    nrow(2)=nrow(2)+1
                    print 104,ok(5:8)
104                 format(16x,4I4)
                    call row3()
                    endif                   
                
                endif
                
            enddo !i3
            endif
          
       enddo ! i2
       endif
    enddo !i1
    end subroutine row2
    subroutine row3()
    integer i1,i2,i3,i4
    do i1=1,upr
        ok(9)=i1
        if(chek(8,i1))then
        
        do i2=1,upr
            ok(10)=i2
            if(chek(9,i2))then
           
            do i3=1,upr                
                ok(11)=i3
                if(chek(10,i3))then
                i4=36-i1-i2-i3                
                    if(chek(11,I4).and. i4>0  .and. i4<17)then                        
                    nrow(3)=nrow(3)+1
                    ok(12)=i4
                    print 104,"three rows:",ok(1:12)
104                 format(a/(20X,4I4))

                    call row4()
                    endif                   
               
                endif
                
            enddo !i3
            endif
          
       enddo ! i2
       endif
    enddo !i1
    end subroutine row3
    subroutine row4()
    integer row4t,s(4)/35,36,37,33/
    i13=s(1)-ok(1)-ok(5)-ok(9) 
    i14=s(2)-ok(2)-ok(6)-ok(10)
    i15=s(3)-ok(3)-ok(7)-ok(11)
    i16=s(4)-ok(4)-ok(8)-ok(12)
    
    if(i13<1 .or. i13>16)go to 90
    if(i14<1 .or. i14>16)go to 90
    if(i15<1 .or. i15>16)go to 90
    if(i16<1 .or. i16>16)go to 90
    print *,"in range:",i13,i14,i15,i16
    read(*,*)
    if(.not.chek(12,i13))go to 90
    ok(13)=i13    
    if(.not.chek(13,i14))go to 90
    ok(14)=i14    
    if(.not.chek(14,i15))go to 90
    ok(15)=i15    
    if(.not.chek(15,i16))go to 90
    ok(16)=i16
    i16a=s(4)-ok(13)-ok(14)-ok(15)
    if(i16.ne.i16a)go to 90
    print 105,"4 X 4 soln:",ok(1:16)
105 format(A/(4I4))
    read(*,*)
90  return
    end subroutine row4
end program