Skip to content

Commit

Permalink
skip ci
Browse files Browse the repository at this point in the history
  • Loading branch information
ctkelley committed Sep 12, 2023
1 parent 83bc165 commit 22e22f6
Show file tree
Hide file tree
Showing 3 changed files with 10 additions and 10 deletions.
2 changes: 1 addition & 1 deletion Codes_For_Docs/Quality.jl
Original file line number Diff line number Diff line change
Expand Up @@ -41,12 +41,12 @@ for idim=1:m
IRData[idim,:]=[N nerra nerre nresa nrese TLP TMP]
# println("errLP = $nerra, errMP=$nerre, ResLP=$nresa, ResMP=$nrese")
end
println(itcount)
if texok
headers = ["N", "ELP", "EMP", "RLP", "RMP", "TLP", "TMP"]
formats = ("%5d & %5.1e & %5.1e & %5.1e & %5.1e & %5.1e & %5.1e")
fprintTeX(headers,formats,IRData)
else
println(itcount)
@printf("%s %s %s %s %s %s %s \n", " N", " ELP",
" EMP", " RLP", " RMP"," TLP", " TMP")
for idim=1:m
Expand Down
4 changes: 2 additions & 2 deletions docs/src/Details.md
Original file line number Diff line number Diff line change
Expand Up @@ -42,7 +42,7 @@ for double, single, and half precision, and the ratio of the half
precision timings to double. The timings came from Julia 1.10-beta2
running on an Apple M2 Pro with 8 performance cores.

Half precision is also difficult to use properly. The low precsion can
Half precision is also difficult to use properly. The low precision can
make iterative refinement fail because the half precision factorization
can have a large error. Here is an example to illustrate this point.
The matrix here is modestly ill-conditioned and you can see that in the
Expand Down Expand Up @@ -78,5 +78,5 @@ julia> norm(z-xd,Inf)
2.34975e-01
```
So you get very poor, but unsurprising, results. While __MultiPrecisionArrays.jl__ supports half precision and I use it all the time, it is not something you would use in your own
work without looking at the literature and making certain you are prepared for strange results. Gettting good results consistently from half precision is an active research area.
work without looking at the literature and making certain you are prepared for strange results. Getting good results consistently from half precision is an active research area.

14 changes: 7 additions & 7 deletions docs/src/index.md
Original file line number Diff line number Diff line change
Expand Up @@ -4,10 +4,10 @@

[MultiPrecisionArrays.jl](https://github.com/ctkelley/MultiPrecisionArrays.jl) is a package for iterative refinement.

This package provides data atructures and solvers for several variants of iterative refinement. It will become much more useful when half precision (aka ```Float16```) is fully supported in LAPACK/BLAS. For now, it's only general-purpose
This package provides data structures and solvers for several variants of iterative refinement. It will become much more useful when half precision (aka ```Float16```) is fully supported in LAPACK/BLAS. For now, it's only general-purpose
application is classical iterative refinement with double precision equations and single precision factorizations.

The half precision stuff is good for those of us doing research in this field. Half precision performace has progressed to the point where you can acutally get things done. On an Apple M2-Pro, a half precision LU only costs 3--5 times
The half precision stuff is good for those of us doing research in this field. Half precision performance has progressed to the point where you can actually get things done. On an Apple M2-Pro, a half precision LU only costs 3--5 times
what a double precision LU costs. This may be as good as it gets unless someone wants to duplicate the LAPACK implementation and get the benefits from blocking, recursion, and clever cache management.

We use a hack-job LU factorization for half precision. Look at the source
Expand All @@ -25,15 +25,15 @@ __IR(A, b)__
- While $\| r \|$ is too large
- Compute the defect $d = (LU)^{-1} r$
- Correct the solution $x = x + d$
- Update the residuial $r = b - Ax$
- Update the residual $r = b - Ax$
- end

In Julia, a code to do this would solve the linear system $A x = b$ in double precision by using a
factorization in a lower precision, say single, within a residual correction iteration. This means that one would need
to allocate storage for a copy of $A$ is the lower precision and factor that copy.

Then one has to determine what the line
$d = (LU)^{-1} r$ means. Do you cast $r$ into the lower precison before the solve or not? __MultiPrecisionArrays.jl__ provides
$d = (LU)^{-1} r$ means. Do you cast $r$ into the lower precision before the solve or not? __MultiPrecisionArrays.jl__ provides
data structures and solvers to manage this.

Here's a simple Julia function for IR that
Expand Down Expand Up @@ -74,15 +74,15 @@ end
The __MPArray__ structure contains both $A$, the low precision copy,
and a vector for the residual.
This lets you allocate the data in advance an reuse the structure
for other right hand sides without rebuilting (or refactoring!) the
for other right hand sides without rebuilding (or refactoring!) the
low precision copy.

As written in the function, the defect uses ```ldiv!``` to compute
```AF\r```. This means that the two triangular factors are stored in
single precision and interprecision transfers are done with each
step in the factorization. While that ```on the fly``` interprecision
transfer is an option, and is needed in many situtations, the
default is to downcase $r$ to low precision, do the solve entirely in
transfer is an option, and is needed in many situations, the
default is to downcast $r$ to low precision, do the solve entirely in
low precision, and the upcast the result. The code for that looks like
```
normr=norm(r)
Expand Down

0 comments on commit 22e22f6

Please sign in to comment.