Cse 245 computer aided circuit simulation and verification
This presentation is the property of its rightful owner.
Sponsored Links
1 / 71

CSE 245: Computer Aided Circuit Simulation and Verification PowerPoint PPT Presentation


  • 84 Views
  • Uploaded on
  • Presentation posted in: General

CSE 245: Computer Aided Circuit Simulation and Verification. Fall 2004, Oct 19 Lecture 7: Matrix Solver II -Iterative Method. Outline. Iterative Method Stationary Iterative Method (SOR, GS,Jacob) Krylov Method (CG, GMRES) Multigrid Method. Iterative Methods. Stationary:

Download Presentation

CSE 245: Computer Aided Circuit Simulation and Verification

An Image/Link below is provided (as is) to download presentation

Download Policy: Content on the Website is provided to you AS IS for your information and personal use and may not be sold / licensed / shared on other websites without getting consent from its author.While downloading, if for some reason you are not able to download a presentation, the publisher may have deleted the file from their server.


- - - - - - - - - - - - - - - - - - - - - - - - - - E N D - - - - - - - - - - - - - - - - - - - - - - - - - -

Presentation Transcript


Cse 245 computer aided circuit simulation and verification

CSE 245: Computer Aided Circuit Simulation and Verification

Fall 2004, Oct 19

Lecture 7:

Matrix Solver II

-Iterative Method


Outline

Outline

  • Iterative Method

    • Stationary Iterative Method (SOR, GS,Jacob)

    • Krylov Method (CG, GMRES)

    • Multigrid Method

Zhengyong (Simon) Zhu, UCSD


Iterative methods

Iterative Methods

Stationary:

x(k+1)=Gx(k)+c

where G and c do not depend on iteration count (k)

Non Stationary:

x(k+1)=x(k)+akp(k)

where computation involves information that change at each iteration

courtesy Alessandra Nardi, UCB


Stationary jacobi method

Stationary: Jacobi Method

In the i-th equation solve for the value of xi while assuming the other entries of x remain fixed:

In matrix terms the method becomes:

where D, -L and -U represent the diagonal, the strictly lower-trg and strictly upper-trg parts of M

M=D-L-U

courtesy Alessandra Nardi, UCB


Stationary gause seidel

Stationary-Gause-Seidel

Like Jacobi, but now assume that previously computed results are used as soon as they are available:

In matrix terms the method becomes:

where D, -L and -U represent the diagonal, the strictly lower-trg and strictly upper-trg parts of M

M=D-L-U

courtesy Alessandra Nardi, UCB


Stationary successive overrelaxation sor

Stationary: Successive Overrelaxation (SOR)

Devised by extrapolation applied to Gauss-Seidel in the form of weighted average:

In matrix terms the method becomes:

where D, -L and -U represent the diagonal, the strictly lower-trg and strictly upper-trg parts of M

M=D-L-U

courtesy Alessandra Nardi, UCB


Cse 245 computer aided circuit simulation and verification

SOR

  • Choose w to accelerate the convergence

    • W =1 : Jacobi / Gauss-Seidel

    • 2>W>1: Over-Relaxation

    • W < 1: Under-Relaxation

Zhengyong (Simon) Zhu, UCSD


Convergence of stationary method

Convergence of Stationary Method

  • Linear Equation: MX=b

  • A sufficient condition for convergence of the solution(GS,Jacob) is that the matrix M is diagonally dominant.

  • If M is symmetric positive definite, SOR converges for any w (0<w<2)

  • A necessary and sufficient condition for the convergence is the magnitude of the largest eigenvalue of the matrix G is smaller than 1

    • Jacobi:

    • Gauss-Seidel

    • SOR:

Zhengyong (Simon) Zhu, UCSD


Outline1

Outline

  • Iterative Method

    • Stationary Iterative Method (SOR, GS,Jacob)

    • Krylov Method (CG, GMRES)

      • Steepest Descent

      • Conjugate Gradient

      • Preconditioning

    • Multigrid Method

Zhengyong (Simon) Zhu, UCSD


Linear equation an optimization problem

Linear Equation: an optimization problem

  • Quadratic function of vector x

  • Matrix A is positive-definite, if for any nonzero vector x

  • If A is symmetric, positive-definite, f(x) is minimized by the solution

Zhengyong (Simon) Zhu, UCSD


Linear equation an optimization problem1

Linear Equation: an optimization problem

  • Quadratic function

  • Derivative

  • If A is symmetric

  • If A is positive-definite

    is minimized by setting to 0

Zhengyong (Simon) Zhu, UCSD


For symmetric positive definite matrix a

For symmetric positive definite matrix A

from J. R. Shewchuk "painless CG"


Gradient of quadratic form

Gradient of quadratic form

The points in the direction of steepest increase of f(x)

from J. R. Shewchuk "painless CG"


Cse 245 computer aided circuit simulation and verification

Symmetric Positive-Definite Matrix A

  • If A is symmetric positive definite

    • P is the arbitrary point

    • X is the solution point

since

We have,

If p != x

Zhengyong (Simon) Zhu, UCSD


If a is not positive definite

If A is not positive definite

  • Positive definite matrix b) negative-definite matrix

  • c) Singular matrix d) positive indefinite matrix

from J. R. Shewchuk "painless CG"


Non stationary iterative method

Non-stationary Iterative Method

  • State from initial guess x0, adjust it until close enough to the exact solution

  • How to choose direction and step size?

i=0,1,2,3,……

Adjustment Direction

Step Size

Zhengyong (Simon) Zhu, UCSD


Steepest descent method 1

Steepest Descent Method (1)

  • Choose the direction in which f decrease most quickly: the direction opposite of

  • Which is also the direction of residue

Zhengyong (Simon) Zhu, UCSD


Steepest descent method 2

Steepest Descent Method (2)

  • How to choose step size ?

    • Line Search

      should minimize f, along the direction of , which means

Orthogonal

Zhengyong (Simon) Zhu, UCSD


Steepest descent algorithm

Steepest Descent Algorithm

Given x0, iterate until residue is smaller than error tolerance

Zhengyong (Simon) Zhu, UCSD


Steepest descent method example

Steepest Descent Method: example

  • Starting at (-2,-2) take the

  • direction of steepest descent of f

  • b) Find the point on the intersec-

  • tion of these two surfaces that

  • minimize f

  • c) Intersection of surfaces.

  • d) The gradient at the bottommost

  • point is orthogonal to the gradient

  • of the previous step

from J. R. Shewchuk "painless CG"


Iterations of steepest descent method

Iterations of Steepest Descent Method

from J. R. Shewchuk "painless CG"


Convergence of steepest descent 1

Convergence of Steepest Descent-1

let

Eigenvector:

EigenValue:

j=1,2,…,n

Energy norm:

Zhengyong (Simon) Zhu, UCSD


Convergence of steepest descent 2

Convergence of Steepest Descent-2

Zhengyong (Simon) Zhu, UCSD


Convergence study n 2

Convergence Study (n=2)

assume

let

Spectral condition number

let

Zhengyong (Simon) Zhu, UCSD


Plot of w

Plot of w

from J. R. Shewchuk "painless CG"


Case study

Case Study

from J. R. Shewchuk "painless CG"


Bound of convergence

Bound of Convergence

It can be proved that it is

also valid for n>2, where

from J. R. Shewchuk "painless CG"


Conjugate gradient method

Conjugate Gradient Method

  • Steepest Descent

    • Repeat search direction

  • Why take exact one step for each direction?

Search direction of Steepest

descent method

figure from J. R. Shewchuk "painless CG"


Orthogonal direction

Orthogonal Direction

Pick orthogonal search direction:

  • We don’t know !!!

Zhengyong (Simon) Zhu, UCSD


Orthogonal a orthogonal

Orthogonal  A-orthogonal

  • Instead of orthogonal search direction, we make search direction A –orthogonal (conjugate)

from J. R. Shewchuk "painless CG"


Search step size

Search Step Size

Zhengyong (Simon) Zhu, UCSD


Iteration finish in n steps

Iteration finish in n steps

Initial error:

A-orthogonal

The error component at direction dj

is eliminated at step j. After n steps,

all errors are eliminated.

Zhengyong (Simon) Zhu, UCSD


Conjugate search direction

Conjugate Search Direction

  • How to construct A-orthogonal search directions, given a set of n linear independent vectors.

  • Since the residue vector in steepest descent method is orthogonal, a good candidate to start with

Zhengyong (Simon) Zhu, UCSD


Construct search direction 1

Construct Search Direction -1

  • In Steepest Descent Method

    • New residue is just a linear combination of previous residue and

  • Let

We have

Krylov SubSpace: repeatedly applying a matrix to a vector

Zhengyong (Simon) Zhu, UCSD


Construct search direction 2

Construct Search Direction -2

let

For i > 0

Zhengyong (Simon) Zhu, UCSD


Construct search direction 3

Construct Search Direction -3

  • can get next direction from the previous one, without saving them all.

let

then

Zhengyong (Simon) Zhu, UCSD


Conjugate gradient algorithm

Conjugate Gradient Algorithm

Given x0, iterate until residue is smaller than error tolerance

Zhengyong (Simon) Zhu, UCSD


Conjugate gradient convergence

Conjugate gradient: Convergence

  • In exact arithmetic, CG converges in n steps (completely unrealistic!!)

  • Accuracy after k steps of CG is related to:

    • consider polynomials of degree k that are equal to 1 at 0.

    • how small can such a polynomial be at all the eigenvalues of A?

  • Thus, eigenvalues close together are good.

  • Condition number:κ(A) = ||A||2 ||A-1||2 = λmax(A) / λmin(A)

  • Residual is reduced by a constant factor by O(κ1/2(A)) iterations of CG.

courtesy J.R.Gilbert, UCSB


Other krylov subspace methods

Other Krylov subspace methods

  • Nonsymmetric linear systems:

    • GMRES: for i = 1, 2, 3, . . . find xi  Ki (A, b) such that ri = (Axi– b)  Ki (A, b)But, no short recurrence => save old vectors => lots more space (Usually “restarted” every k iterations to use less space.)

    • BiCGStab, QMR, etc.:Two spaces Ki (A, b)and Ki (AT, b)w/ mutually orthogonal bases Short recurrences => O(n) space, but less robust

    • Convergence and preconditioning more delicate than CG

    • Active area of current research

  • Eigenvalues: Lanczos (symmetric), Arnoldi (nonsymmetric)

courtesy J.R.Gilbert, UCSB


Preconditioners

Preconditioners

  • Suppose you had a matrix B such that:

    • condition number κ(B-1A) is small

    • By = z is easy to solve

  • Then you could solve (B-1A)x = B-1b instead of Ax = b

  • B = A is great for (1), not for (2)

  • B = I is great for (2), not for (1)

  • Domain-specific approximations sometimes work

  • B = diagonal of A sometimes works

  • Better: blend in some direct-methods ideas. . .

courtesy J.R.Gilbert, UCSB


Preconditioned conjugate gradient iteration

Preconditioned conjugate gradient iteration

  • One matrix-vector multiplication per iteration

  • One solve with preconditioner per iteration

x0 = 0, r0 = b, d0 = B-1r0, y0 = B-1r0

for k = 1, 2, 3, . . .

αk = (yTk-1rk-1) / (dTk-1Adk-1) step length

xk = xk-1 + αk dk-1 approx solution

rk = rk-1 – αk Adk-1 residual

yk = B-1rk preconditioning solve

βk = (yTk rk) / (yTk-1rk-1) improvement

dk = yk + βk dk-1 search direction

courtesy J.R.Gilbert, UCSB


Outline2

Outline

  • Iterative Method

    • Stationary Iterative Method (SOR, GS,Jacob)

    • Krylov Method (CG, GMRES)

    • Multigrid Method

Zhengyong (Simon) Zhu, UCSD


What is the multigrid

What is the multigrid

  • A multilevel iterative method to solve

    • Ax=b

  • Originated in PDEs on geometric grids

  • Expend the multigrid idea to unstructured problem – Algebraic MG

  • Geometric multigrid for presenting the basic ideas of the multigrid method.

Zhengyong (Simon) Zhu, UCSD


The model problem

v3

v4

v1

v2

v5

v6

v8

v7

+

vs

The model problem

Ax = b

Zhengyong (Simon) Zhu, UCSD


Simple iterative method

Simple iterative method

  • x(0) -> x(1) -> … -> x(k)

  • Jacobi iteration

  • Matrix form : x(k) = Rjx(k-1) + Cj

  • General form: x(k) = Rx(k-1) + C (1)

  • Stationary: x* = Rx* + C (2)

Zhengyong (Simon) Zhu, UCSD


Error and convergence

Error and Convergence

Definition: errore = x* - x (3)

residualr = b – Ax (4)

e, r relation: Ae = r (5) ((3)+(4))

e(1) = x*-x(1) = Rx* + C – Rx(0)– C =Re(0)

Error equatione(k) = Rke(0) (6) ((1)+(2)+(3))

Convergence:

Zhengyong (Simon) Zhu, UCSD


Error of diffenent frequency

k= 1

k= 4

k= 2

Error of diffenent frequency

  • Wavenumber k and frequency 

  • = k/n

  • High frequency error is more oscillatory between points

Zhengyong (Simon) Zhu, UCSD


Iteration reduce low frequency error efficiently

Iteration reduce low frequency error efficiently

  • Smoothing iteration reduce high frequency error efficiently, but not low frequency error

Error

k = 1

k = 2

k = 4

Iterations

Zhengyong (Simon) Zhu, UCSD


Multigrid a first glance

2

1

3

4

3

4

1

2

5

6

8

7

Multigrid – a first glance

  • Two levels : coarse and fine grid

2h

A2hx2h=b2h

h

Ahxh=bh

Ax=b

Zhengyong (Simon) Zhu, UCSD


Idea 1 the v cycle iteration

Idea 1: the V-cycle iteration

  • Also called the nested iteration

Start with

2h

A2hx2h = b2h

A2hx2h = b2h

Iterate =>

Prolongation: 

Restriction: 

h

Ahxh = bh

Iterate to get

Question 1: Why we need the coarse grid ?

Zhengyong (Simon) Zhu, UCSD


Prolongation

2

1

3

4

3

4

1

2

5

6

8

7

Prolongation

  • Prolongation (interpolation) operator

    xh = x2h

Zhengyong (Simon) Zhu, UCSD


Restriction

2

1

3

4

3

4

1

2

5

6

8

7

Restriction

  • Restriction operator

    xh = x2h

Zhengyong (Simon) Zhu, UCSD


Smoothing

Smoothing

  • The basic iterations in each level

    In ph: xphold  xphnew

  • Iteration reduces the error, makes the error smooth geometrically.

    So the iteration is called smoothing.

Zhengyong (Simon) Zhu, UCSD


Why multilevel

Why multilevel ?

  • Coarse lever iteration is cheap.

  • More than this…

    • Coarse level smoothing reduces the error more efficiently than fine level in some way .

    • Why ? ( Question 2 )

Zhengyong (Simon) Zhu, UCSD


Error restriction

Error restriction

  • Map error to coarse grid will make the error more oscillatory

K = 4,  = 

K = 4,  = /2

Zhengyong (Simon) Zhu, UCSD


Idea 2 residual correction

Idea 2: Residual correction

  • Known current solution x

  • Solve Ax=b eq. to

  • MG do NOT map x directly between levels

    Map residual equation to coarse level

    • Calculate rh

    • b2h= Ih2h rh ( Restriction )

    • eh =Ih2hx2h ( Prolongation )

    • xh = xh + eh

Zhengyong (Simon) Zhu, UCSD


Why residual correction

Why residual correction ?

  • Error is smooth at fine level, but the actual solution may not be.

  • Prolongation results in a smooth error in fine level, which is suppose to be a good evaluation of the fine level error.

  • If the solution is not smooth in fine level, prolongation will introduce more high frequency error.

Zhengyong (Simon) Zhu, UCSD


Revised v cycle with idea 2

Revised V-cycle with idea 2

  • Smoothing on xh

  • Calculate rh

  • b2h= Ih2h rh

  • Smoothing on x2h

  • eh =Ih2hx2h

  • Correct: xh = xh + eh

2h h

`

Restriction

Prolongation

Zhengyong (Simon) Zhu, UCSD


What is a 2h

What is A2h

  • Galerkin condition

Zhengyong (Simon) Zhu, UCSD


Going to multilevels

Going to multilevels

  • V-cycle and W-cycle

  • Full Multigrid V-cycle

h

2h

4h

h

2h

4h

8h

Zhengyong (Simon) Zhu, UCSD


Performance of multigrid

Performance of Multigrid

  • Complexity comparison

Zhengyong (Simon) Zhu, UCSD


Summary of mg ideas

Summary of MG ideas

Three important ideas of MG

  • Nested iteration

  • Residual correction

  • Elimination of error:

    high frequency : fine grid

    low frequency : coarse grid

Zhengyong (Simon) Zhu, UCSD


Amg for unstructured grids

1

2

4

3

5

6

AMG :for unstructured grids

  • Ax=b, no regular grid structure

  • Fine grid defined from A

Zhengyong (Simon) Zhu, UCSD


Three questions for amg

Three questions for AMG

  • How to choose coarse grid

  • How to define the smoothness of errors

  • How are interpolation and prolongation done

Zhengyong (Simon) Zhu, UCSD


How to choose coarse grid

How to choose coarse grid

  • Idea:

    • C/F splitting

    • As few coarse grid point as possible

    • For each F-node, at least one of its neighbor is a C-node

    • Choose node with strong coupling to other nodes as C-node

1

2

4

3

5

6

Zhengyong (Simon) Zhu, UCSD


How to define the smoothness of error

How to define the smoothness of error

  • AMG fundamental concept:

    Smooth error = small residuals

  • ||r|| << ||e||

Zhengyong (Simon) Zhu, UCSD


How are prolongation and restriction done

How are Prolongation and Restriction done

  • Prolongation is based on smooth error and strong connections

  • Common practice: I

Zhengyong (Simon) Zhu, UCSD


Amg prolongation 2

AMG Prolongation (2)

Zhengyong (Simon) Zhu, UCSD


Amg prolongation 3

AMG Prolongation (3)

  • Restriction :

Zhengyong (Simon) Zhu, UCSD


Summary

Summary

  • Multigrid is a multilevel iterative method.

  • Advantage: scalable

  • If no geometrical grid is available, try Algebraic multigrid method

Zhengyong (Simon) Zhu, UCSD


The landscape of solvers

Direct

A = LU

Iterative

y’ = Ay

More General

Non-

symmetric

Symmetric

positive

definite

More Robust

The landscape of Solvers

More Robust

Less Storage (if sparse)

courtesy J.R.Gilbert, UCSB


  • Login