Skip to content

Commit

Permalink
[bug] Fix MatrixFreeCG so it can handle multiple input sizes. (#8070)
Browse files Browse the repository at this point in the history
### Brief Summary

Previously, `MatrixFreeCG` inappropriately assume that the intermediate
vectors used in CG solver are all two-dimensional `ti.field`. This is
incorrect because the dimension of `Ap` `p` and `r` should be consistent
with input `b` and `x`.

### Walkthrough
Instead of simply put `vector_fields_builder.dense(ti.ij, size).place(p,
r, Ap)`, now it's based on the length of `size`.
```python
def MatrixFreeCG(A, b, x, tol=1e-6, maxiter=5000, quiet=True):
    ...                                                                                                                                                                          
    vector_fields_builder = ti.FieldsBuilder()                                                                                                                                                                                          
    p = ti.field(dtype=solver_dtype)                                                                                                                                                                                                    
    r = ti.field(dtype=solver_dtype)                                                                                                                                                                                                    
    Ap = ti.field(dtype=solver_dtype)                                                                                                                                                                                                   
    if len(size) == 1:     # Determine the `axes` argument based on the length of `size`                                                                                                                                                                                                             
        axes = ti.i                                                                                                                                                                                                                     
    elif len(size) == 2:                                                                                                                                                                                                                
        axes = ti.ij                                                                                                                                                                                                                    
    elif len(size) == 3:                                                                                                                                                                                                                
        axes = ti.ijk                                                                                                                                                                                                                   
    else:                                                                                                                                                                                                                               
        raise TaichiRuntimeError(f"MatrixFreeCG currently cannot support {len(size)}-D inputs.")                                                                                                                                        
    vector_fields_builder.dense(axes, size).place(p, r, Ap) 
```
  • Loading branch information
houkensjtu authored May 31, 2023
1 parent e3cccb8 commit 3332eee
Showing 1 changed file with 9 additions and 1 deletion.
10 changes: 9 additions & 1 deletion python/taichi/linalg/matrixfree_cg.py
Original file line number Diff line number Diff line change
Expand Up @@ -47,7 +47,15 @@ def MatrixFreeCG(A, b, x, tol=1e-6, maxiter=5000, quiet=True):
p = ti.field(dtype=solver_dtype)
r = ti.field(dtype=solver_dtype)
Ap = ti.field(dtype=solver_dtype)
vector_fields_builder.dense(ti.ij, size).place(p, r, Ap)
if len(size) == 1:
axes = ti.i
elif len(size) == 2:
axes = ti.ij
elif len(size) == 3:
axes = ti.ijk
else:
raise TaichiRuntimeError(f"MatrixFreeCG only support 1D, 2D, 3D inputs; your inputs is {len(size)}-D.")
vector_fields_builder.dense(axes, size).place(p, r, Ap)
vector_fields_snode_tree = vector_fields_builder.finalize()

scalar_builder = ti.FieldsBuilder()
Expand Down

0 comments on commit 3332eee

Please sign in to comment.