anchor.diff:     u=z(:,j); v=z(ip,:)-z(i,:)
      zinv=(1.d0,0.d0)/z(i,j)
   
  -!$OMP parallel do shared(z,zinv,u,v)
      do jj=1,n
       call zaxpy(m,zinv*v(jj),u,1,z(1,jj),1)
      end do
  -!$OMP END parallel do
      z(i,:)=(0.d0,0.d0); z(i,j)=(1.d0,0.d0)   ! FORCE
     end do 
   100 continue
  
positive.code: 
  !!  call dger(m,n,zinv,u,1,v,1,z,m)
  
  !!..OMP PARALLEL WORKSHARE SHARED(z,u,v,zinv) 
  !!   forall(ii=1:m,jj=1:n) z(ii,jj)=z(ii,jj)+zinv*u(ii)*v(jj)
  !!..OMP END PARALLEL WORKSHARE
  
  !$OMP parallel do shared(z,zinv,u,v)
     do jj=1,n
      call daxpy(m,zinv*v(jj),u,1,z(1,jj),1)
     end do
  !$OMP END parallel do
     z(i,:)=0.d0; z(i,j)=1.d0   ! FORCE
    end do 
  100 continue
  
negative.code: !   allocate(tt(m,n)); call maxvola_d(a,p,tt); nrm=dnrm2(m*n,tt,1); tt=tt-z; err=dnrm2(m*n,tt,1); deallocate(tt)
  !   write(*,'(a,i4,3(a,e9.3))')'iter: ',iter,' err: ',err,' nrm: ',nrm,' relerr: ',err/nrm
  !--
  
     do jj=1,n
      q(jj)=idamax(m,z(1,jj),1)
      v(jj)=abs(z(q(jj),jj))
     end do
     j=maxloc(v,1); i=q(j)
     
     amax=dabs(z(i,j))
  