The paper presents an algebraic scheme to construct hierarchical divergence-free basis for velocity in incompressible fluids. A reduced system of equations is solved in the corresponding subspace by an appropriate iterative method. The basis is constructed from the matrix representing the incompressibility constraints by computing algebraic decompositions of local constraint matrices. A recursive strategy leads to a hierarchical basis with desirable properties such as fast matrix-vector products, a well-conditioned reduced system, and efficient parallelization of the computation. The scheme has been extended to particulate flow problems in which the Navier-Stokes equations for fluid are coupled with equations of motion for rigid particles suspended in the fluid. Experimental results of particulate flow simulations have been reported for the SGI Origin 2000.