**FAME-matlab Package: Fast Algorithm for ** **Maxwell Equations **

**Tsung-Ming Huang
**

### Modelling, Simulation and Analysis of

### Nonlinear Optics, NUK, September, 4-8, 2015

### 1

### 2

### 3

**FAME** group

### 2

**Wen-Wei Lin **

### Department of Applied Mathematics National Chiao-Tung University

**Weichung Wang **

### Department of Mathematics National Taiwan University

**Han-En Hsieh(謝函恩) **

### Department of Mathematics National Taiwan University

**Chien-Chih Huang(黃建智) **

### Department of Mathematics

### National Taiwan Normal University

**F** ast **A** lgorithm for **M** axwell **E** q

### 4

### FAME

### FAME. ^{m}

### FAME. ^{mpi}

### FAME. ^{GPU}

### FAME

### 5

### Eigen-solvers (JD, SIRA)

**Dispersive Metallic **

**materials** **Complex materials** **Photonic Crystals**

**Simple Cubic ** **Face-Centered Cubic **

**∇ × ∇ × E(r) =** λε ^{(r)E(r)} **∇ × ∇ × E(r) =** ω ^{2} ε ^{(r,} ω ^{)E(r)} ^{⎡} _{⎣} ^{⎢} ^{∇ ×} _{0} _{∇ ×} ^{0} ^{⎤} _{⎦} ^{⎥} ^{⎡} _{⎣} ^{⎢} _{H} ^{E} ^{⎤} _{⎦} ^{⎥ = i} ^{ω} ^{⎡} _{−} ^{ζ} _{ε} _{−} ^{µ} _{ξ}

**∇ × ∇ × E(r) =**

^{(r)E(r)}**∇ × ∇ × E(r) =**

^{(r,}

^{)E(r)}_{H}

^{E}

### ⎣ ⎢

### ⎢

### ⎤

### ⎦ ⎥

### ⎥ *E* *H*

### ⎡

### ⎣ ⎢ ⎤

### ⎦ ⎥

employing projections onto the range space of the discrete
matrix. (However, it is noted that the problem of many zeros
does not occur in the PWE as in each plane wave component
the number of polarizations is chosen to be the independent
two rather three.) (iii) As the third issue of importance, it is
of great interest to develop a fast algorithm for the
eigenvalue problem. An inverse method, accelerated by
multigrid technique with use of projection is proposed for
this purpose. The method exploits the sparsity of the matrix
for the eigenvalue problem in the finite diﬀerence formu-
lation. Because of the above mentioned diﬃculties, the
present method is a nontrivial extension of a similar method,
recently developed by the authors^{24)}for computing photonic
band structures in two dimensions.

In this study, we compute the band structures for three
types of photonic structures. The first one is a modified
simple cubic lattice consisting of dielectric spheres on the
lattice sites, each connected to its six nearest neighbors by
thin circular cylinders, which was proposed by Biswas et
al.^{25)}Figure 1 shows the modified simple cubic lattice. It is
noted that the original simple cubic structure comprising a
lattice of rods has been fabricated recently with advanced
silicon processing techniques.^{26)}The second one is the
tetragonal square spiral structure comprising a lattice of
circular or square cylinders, which was proposed by Toader
and John,^{27,28)}as shown in Fig. 2. Spiral structure was
discussed previously by Chutinan and Noda.^{29)}The square
spiral structure is arranged to connect the lattice points of
diamond structure with specific order, and is amenable to the
current technique of fabrication GLAD (GLancing Angle
Deposition) as discussed in refs. 30 and 31. As a third
example, we propose a diamond structure that has sp^{3}-like
configuration, composed of dielectric spheres with connect-
ing spheroids, as shown in Fig. 3. Diamond structures are
known to have large band gaps between relatively lower
branches either in diamond network or inverse diamond
structure.^{12,28)}In the present study, the spheroids, instead of
circular cylinders, take the positions of ‘‘valence bonds’’ to
imitate the sp^{3}structure of the electrons of diamond atoms.

Recently, submicron diamond-lattice photonic crystals have been successfully produced by two-photon laser nanofabri-

cation (photopolymerization).^{32)}

The order of presentation of the paper is organized as
follows. In §2, we show how to correctly formulate the finite
diﬀerence method for the double curl operator of the
photonic eigenvalue problem. In §3, we develop the
numerical method (inverse iteration with the full multigrid
acceleration) and present the fast algorithm, in which two
alternative methods of projection are proposed to avoid the
necessity of deflating zeros). In §4, we first present
numerical results that illustrate the eﬃciency of the
presently developed method. Then, the band structures are
computed for the modified simple cubic lattice, the tetrag-
onal square spiral structure (direct and inverse structure) and
the diamond structure with sp^{3}-like configuration. Finally,
concluding remarks with a summary of results are drawn in

§5.

2. Basic Equations and Finite Diﬀerence Formulation The electromagnetic waves propagating in the photonic crystals are well described by Maxwell’s equations. For linear isotropic and frequency-independent dielectric mate- rials with permeability close to one, the time-harmonic Fig. 1. Modified simple cubic lattice comprising dielectric spheres and

connecting thin circular cylinders.^{25)}

Fig. 2. Tetragonal square spiral structure comprising circular cylin-
ders.^{27,28)}

Fig. 3. diamond structure with sp^{3}-like configuration comprising dielec-
tric spheres and connecting spheroids.

728 J. Phys. Soc. Jpn., Vol. 73, No. 3, March, 2004 R. L. CHERNet al.

employing projections onto the range space of the discrete
matrix. (However, it is noted that the problem of many zeros
does not occur in the PWE as in each plane wave component
the number of polarizations is chosen to be the independent
two rather three.) (iii) As the third issue of importance, it is
of great interest to develop a fast algorithm for the
eigenvalue problem. An inverse method, accelerated by
multigrid technique with use of projection is proposed for
this purpose. The method exploits the sparsity of the matrix
for the eigenvalue problem in the finite diﬀerence formu-
lation. Because of the above mentioned diﬃculties, the
present method is a nontrivial extension of a similar method,
recently developed by the authors^{24)}for computing photonic
band structures in two dimensions.

In this study, we compute the band structures for three
types of photonic structures. The first one is a modified
simple cubic lattice consisting of dielectric spheres on the
lattice sites, each connected to its six nearest neighbors by
thin circular cylinders, which was proposed by Biswas et
al.^{25)}Figure 1 shows the modified simple cubic lattice. It is
noted that the original simple cubic structure comprising a
lattice of rods has been fabricated recently with advanced
silicon processing techniques.^{26)} The second one is the
tetragonal square spiral structure comprising a lattice of
circular or square cylinders, which was proposed by Toader
and John,^{27,28)} as shown in Fig. 2. Spiral structure was
discussed previously by Chutinan and Noda.^{29)}The square
spiral structure is arranged to connect the lattice points of
diamond structure with specific order, and is amenable to the
current technique of fabrication GLAD (GLancing Angle
Deposition) as discussed in refs. 30 and 31. As a third
example, we propose a diamond structure that has sp^{3}-like
configuration, composed of dielectric spheres with connect-
ing spheroids, as shown in Fig. 3. Diamond structures are
known to have large band gaps between relatively lower
branches either in diamond network or inverse diamond
structure.^{12,28)}In the present study, the spheroids, instead of
circular cylinders, take the positions of ‘‘valence bonds’’ to
imitate the sp^{3}structure of the electrons of diamond atoms.

Recently, submicron diamond-lattice photonic crystals have been successfully produced by two-photon laser nanofabri-

cation (photopolymerization).^{32)}

The order of presentation of the paper is organized as
follows. In §2, we show how to correctly formulate the finite
diﬀerence method for the double curl operator of the
photonic eigenvalue problem. In §3, we develop the
numerical method (inverse iteration with the full multigrid
acceleration) and present the fast algorithm, in which two
alternative methods of projection are proposed to avoid the
necessity of deflating zeros). In §4, we first present
numerical results that illustrate the eﬃciency of the
presently developed method. Then, the band structures are
computed for the modified simple cubic lattice, the tetrag-
onal square spiral structure (direct and inverse structure) and
the diamond structure with sp^{3}-like configuration. Finally,
concluding remarks with a summary of results are drawn in

§5.

2. Basic Equations and Finite Diﬀerence Formulation The electromagnetic waves propagating in the photonic crystals are well described by Maxwell’s equations. For linear isotropic and frequency-independent dielectric mate- rials with permeability close to one, the time-harmonic Fig. 1. Modified simple cubic lattice comprising dielectric spheres and

connecting thin circular cylinders.^{25)}

Fig. 2. Tetragonal square spiral structure comprising circular cylin-
ders.^{27,28)}

Fig. 3. diamond structure with sp^{3}-like configuration comprising dielec-
tric spheres and connecting spheroids.

728 J. Phys. Soc. Jpn., Vol. 73, No. 3, March, 2004 R. L. CHERNet al.

**Generalized eigenvalue ** **problems for 3D **

**photonic crystal**

### 6

### Curl operator

### Central edge points

### Central face points

### where

### Resulting generalized eigenvalue problem

### with diagonal B

### 7

**∇ × ∇ × E(r) =** ω ^{2} ε ^{(r)E(r)}

**∇ × ∇ × E(r) =**

^{(r)E(r)}**∇ × E(r) = H(r) ⇒ Ce = h**

**∇ × E(r) = H(r) ⇒ Ce = h**

**∇ × H(r) =** ω ^{2} ε **(r)E(r)** *⇒ C* ^{∗} **h** = ω ^{2} **B e**

**∇ × H(r) =**

**(r)E(r)**

**B e**

*∇ × E =*

### 0 − ∂

*∂z*

### ∂

*∂y*

### ∂

*∂z* 0 − ∂

*∂x*

### − ∂ *∂y*

### ∂

*∂x* 0

### ⎡

### ⎣

### ⎢ ⎢

### ⎢ ⎢

### ⎢ ⎢

### ⎢

### ⎤

### ⎦

### ⎥ ⎥

### ⎥ ⎥

### ⎥ ⎥

### ⎥ *E*

_{1}

*E*

_{2}

*E*

_{3}

### ⎡

### ⎣

### ⎢ ⎢

### ⎢

### ⎤

### ⎦

### ⎥ ⎥

### ⎥

*C* =

### 0 *−C* _{3} *C* _{2} *C* _{3} 0 *−C* _{1}

*−C* _{2} *C* _{1} 0

### ⎡

### ⎣

### ⎢ ⎢

### ⎢

### ⎤

### ⎦

### ⎥ ⎥

### ⎥ ∈! ^{3n×3n}

^{3n×3n}

*C* _{1} *= I* _{n} _{2} _{n} _{3} *⊗ K* _{1} ∈! ^{n} ^{×n} *, C* _{2} *= I* _{n} _{3} *⊗ K* _{2} ∈! ^{n} ^{×n} *, C* _{3} *= K* _{3} ∈! ^{n} ^{×n}

_{n}

_{n}

^{n}

^{×n}

_{n}

^{n}

^{×n}

^{n}

^{×n}

*C* ^{∗} *C* − ω ^{2} *B*

### ( ) ^{x} ^{≡ A −} ^{(} ^{λ} ^{B} ^{)} ^{x} ^{= 0}

^{x}

^{≡ A −}

^{B}

^{x}

**Finite Diff. Assoc. with Quasi-Periodic Cond.**

### 8

*K* _{1} = 1 δ _{x}

_{x}

### −1 1

### ! ! −1 1 *e* ^{ı2} ^{π} ^{k⋅a} ^{1} −1

^{ı2}

^{k⋅a}

### ⎡

### ⎣

### ⎢ ⎢

### ⎢ ⎢

### ⎤

### ⎦

### ⎥ ⎥

### ⎥ ⎥

### ∈! ^{n} ^{1} ^{×n} ^{1} ,

^{n}

^{×n}

*K* _{2} = 1 δ _{y}

_{y}

*−I* _{n} _{1} *I* _{n} _{1}

_{n}

_{n}

### ! !

*−I* _{n} _{1} *I* _{n}

_{n}

_{n}

### 1

*e* ^{ı2} ^{π} ^{k⋅a} ^{2} *J* _{2} *−I* _{n} _{1}

^{ı2}

^{k⋅a}

_{n}

### ⎡

### ⎣

### ⎢ ⎢

### ⎢ ⎢

### ⎢

### ⎤

### ⎦

### ⎥ ⎥

### ⎥ ⎥

### ⎥

### ∈! ^{(n} ^{1} ^{n} ^{2} ^{)×(n} ^{1} ^{n} ^{2} ^{)} ,

^{(n}

^{n}

^{n}

*K* _{3} = 1 δ _{z}

_{z}

*−I* _{n} _{1} _{n} _{2} *I* _{n} _{1} _{n} _{2}

_{n}

_{n}

_{n}

_{n}

### ! !

*−I* _{n} _{1} _{n} _{2} *I* _{n}

_{n}

_{n}

_{n}

### 1 *n* _{2}

*e* ^{ı2} ^{π} ^{k⋅a} ^{3} *J* _{3} *−I* _{n} _{1} _{n} _{2}

^{ı2}

^{k⋅a}

_{n}

_{n}

### ⎡

### ⎣

### ⎢ ⎢

### ⎢ ⎢

### ⎢

### ⎤

### ⎦

### ⎥ ⎥

### ⎥ ⎥

### ⎥

### ∈! ^{n×n}

^{n×n}

**Finite Diff. Assoc. with Quasi-Periodic Cond.**

### 8

*K* _{1} = 1 δ _{x}

_{x}

### −1 1

### ! ! −1 1 *e* ^{ı2} ^{π} ^{k⋅a} ^{1} −1

^{ı2}

^{k⋅a}

### ⎡

### ⎣

### ⎢ ⎢

### ⎢ ⎢

### ⎤

### ⎦

### ⎥ ⎥

### ⎥ ⎥

### ∈! ^{n} ^{1} ^{×n} ^{1} ,

^{n}

^{×n}

*K* _{2} = 1 δ _{y}

_{y}

*−I* _{n} _{1} *I* _{n} _{1}

_{n}

_{n}

### ! !

*−I* _{n} _{1} *I* _{n}

_{n}

_{n}

### 1

*e* ^{ı2} ^{π} ^{k⋅a} ^{2} *J* _{2} *−I* _{n} _{1}

^{ı2}

^{k⋅a}

_{n}

### ⎡

### ⎣

### ⎢ ⎢

### ⎢ ⎢

### ⎢

### ⎤

### ⎦

### ⎥ ⎥

### ⎥ ⎥

### ⎥

### ∈! ^{(n} ^{1} ^{n} ^{2} ^{)×(n} ^{1} ^{n} ^{2} ^{)} ,

^{(n}

^{n}

^{n}

*K* _{3} = 1 δ _{z}

_{z}

*−I* _{n} _{1} _{n} _{2} *I* _{n} _{1} _{n} _{2}

_{n}

_{n}

_{n}

_{n}

### ! !

*−I* _{n} _{1} _{n} _{2} *I* _{n}

_{n}

_{n}

_{n}

### 1 *n* _{2}

*e* ^{ı2} ^{π} ^{k⋅a} ^{3} *J* _{3} *−I* _{n} _{1} _{n} _{2}

^{ı2}

^{k⋅a}

_{n}

_{n}

### ⎡

### ⎣

### ⎢ ⎢

### ⎢ ⎢

### ⎢

### ⎤

### ⎦

### ⎥ ⎥

### ⎥ ⎥

### ⎥

### ∈! ^{n×n}

^{n×n}

**E(r** **+ a** _{ℓ} ) *= e* ^{i2} ^{π} ^{k⋅a} ^{ℓ} **E(r)**

**E(r**

^{k⋅a}

**E(r)**

**Finite Diff. Assoc. with Quasi-Periodic Cond.**

### 8

*K* _{1} = 1 δ _{x}

_{x}

### −1 1

### ! ! −1 1 *e* ^{ı2} ^{π} ^{k⋅a} ^{1} −1

^{ı2}

^{k⋅a}

### ⎡

### ⎣

### ⎢ ⎢

### ⎢ ⎢

### ⎤

### ⎦

### ⎥ ⎥

### ⎥ ⎥

### ∈! ^{n} ^{1} ^{×n} ^{1} ,

^{n}

^{×n}

*K* _{2} = 1 δ _{y}

_{y}

*−I* _{n} _{1} *I* _{n} _{1}

_{n}

_{n}

### ! !

*−I* _{n} _{1} *I* _{n}

_{n}

_{n}

### 1

*e* ^{ı2} ^{π} ^{k⋅a} ^{2} *J* _{2} *−I* _{n} _{1}

^{ı2}

^{k⋅a}

_{n}

### ⎡

### ⎣

### ⎢ ⎢

### ⎢ ⎢

### ⎢

### ⎤

### ⎦

### ⎥ ⎥

### ⎥ ⎥

### ⎥

### ∈! ^{(n} ^{1} ^{n} ^{2} ^{)×(n} ^{1} ^{n} ^{2} ^{)} ,

^{(n}

^{n}

^{n}

*K* _{3} = 1 δ _{z}

_{z}

*−I* _{n} _{1} _{n} _{2} *I* _{n} _{1} _{n} _{2}

_{n}

_{n}

_{n}

_{n}

### ! !

*−I* _{n} _{1} _{n} _{2} *I* _{n}

_{n}

_{n}

_{n}

### 1 *n* _{2}

*e* ^{ı2} ^{π} ^{k⋅a} ^{3} *J* _{3} *−I* _{n} _{1} _{n} _{2}

^{ı2}

^{k⋅a}

_{n}

_{n}

### ⎡

### ⎣

### ⎢ ⎢

### ⎢ ⎢

### ⎢

### ⎤

### ⎦

### ⎥ ⎥

### ⎥ ⎥

### ⎥

### ∈! ^{n×n}

^{n×n}

**Finite Diff. Assoc. with Quasi-Periodic Cond.**

### 8

*K* _{1} = 1 δ _{x}

_{x}

### −1 1

### ! ! −1 1 *e* ^{ı2} ^{π} ^{k⋅a} ^{1} −1

^{ı2}

^{k⋅a}

### ⎡

### ⎣

### ⎢ ⎢

### ⎢ ⎢

### ⎤

### ⎦

### ⎥ ⎥

### ⎥ ⎥

### ∈! ^{n} ^{1} ^{×n} ^{1} ,

^{n}

^{×n}

*K* _{2} = 1 δ _{y}

_{y}

*−I* _{n} _{1} *I* _{n} _{1}

_{n}

_{n}

### ! !

*−I* _{n} _{1} *I* _{n}

_{n}

_{n}

### 1

*e* ^{ı2} ^{π} ^{k⋅a} ^{2} *J* _{2} *−I* _{n} _{1}

^{ı2}

^{k⋅a}

_{n}

### ⎡

### ⎣

### ⎢ ⎢

### ⎢ ⎢

### ⎢

### ⎤

### ⎦

### ⎥ ⎥

### ⎥ ⎥

### ⎥

### ∈! ^{(n} ^{1} ^{n} ^{2} ^{)×(n} ^{1} ^{n} ^{2} ^{)} ,

^{(n}

^{n}

^{n}

*K* _{3} = 1 δ _{z}

_{z}

*−I* _{n} _{1} _{n} _{2} *I* _{n} _{1} _{n} _{2}

_{n}

_{n}

_{n}

_{n}

### ! !

*−I* _{n} _{1} _{n} _{2} *I* _{n}

_{n}

_{n}

_{n}

### 1 *n* _{2}

*e* ^{ı2} ^{π} ^{k⋅a} ^{3} *J* _{3} *−I* _{n} _{1} _{n} _{2}

^{ı2}

^{k⋅a}

_{n}

_{n}

### ⎡

### ⎣

### ⎢ ⎢

### ⎢ ⎢

### ⎢

### ⎤

### ⎦

### ⎥ ⎥

### ⎥ ⎥

### ⎥

### ∈! ^{n×n}

^{n×n}

*J* _{2} *= I* _{n} _{1} *, J* _{3} *= I* _{n} _{1} _{n} _{2}

_{n}

_{n}

_{n}

### For SC lattice

**Finite Diff. Assoc. with Quasi-Periodic Cond.**

### 8

*K* _{1} = 1 δ _{x}

_{x}

### −1 1

### ! ! −1 1 *e* ^{ı2} ^{π} ^{k⋅a} ^{1} −1

^{ı2}

^{k⋅a}

### ⎡

### ⎣

### ⎢ ⎢

### ⎢ ⎢

### ⎤

### ⎦

### ⎥ ⎥

### ⎥ ⎥

### ∈! ^{n} ^{1} ^{×n} ^{1} ,

^{n}

^{×n}

*K* _{2} = 1 δ _{y}

_{y}

*−I* _{n} _{1} *I* _{n} _{1}

_{n}

_{n}

### ! !

*−I* _{n} _{1} *I* _{n}

_{n}

_{n}

### 1

*e* ^{ı2} ^{π} ^{k⋅a} ^{2} *J* _{2} *−I* _{n} _{1}

^{ı2}

^{k⋅a}

_{n}

### ⎡

### ⎣

### ⎢ ⎢

### ⎢ ⎢

### ⎢

### ⎤

### ⎦

### ⎥ ⎥

### ⎥ ⎥

### ⎥

### ∈! ^{(n} ^{1} ^{n} ^{2} ^{)×(n} ^{1} ^{n} ^{2} ^{)} ,

^{(n}

^{n}

^{n}

*K* _{3} = 1 δ _{z}

_{z}

*−I* _{n} _{1} _{n} _{2} *I* _{n} _{1} _{n} _{2}

_{n}

_{n}

_{n}

_{n}

### ! !

*−I* _{n} _{1} _{n} _{2} *I* _{n}

_{n}

_{n}

_{n}

### 1 *n* _{2}

*e* ^{ı2} ^{π} ^{k⋅a} ^{3} *J* _{3} *−I* _{n} _{1} _{n} _{2}

^{ı2}

^{k⋅a}

_{n}

_{n}

### ⎡

### ⎣

### ⎢ ⎢

### ⎢ ⎢

### ⎢

### ⎤

### ⎦

### ⎥ ⎥

### ⎥ ⎥

### ⎥

### ∈! ^{n×n}

^{n×n}

*J* _{2} *= I* _{n} _{1} *, J* _{3} *= I* _{n} _{1} _{n} _{2}

_{n}

_{n}

_{n}

### For SC lattice

*J* _{2} = 0 *e* ^{−ı2} ^{πk⋅a} ^{1} *I* _{n}

^{−ı2}

^{πk⋅a}

_{n}

### 1 /2

*I* _{n}

_{n}

### 1 /2 0

### ⎡

### ⎣

### ⎢ ⎢

### ⎤

### ⎦

### ⎥ ⎥ ∈! ^{n} ^{1} ^{×n} ^{1} ,

^{n}

^{×n}

*J* _{3} =

### 0 *e* ^{−ı2πk⋅a} ^{2} *I* _{1}

^{−ı2πk⋅a}

### 3 *n* _{2} *⊗ I* _{n}

_{n}

### 1

*I* _{2}

### 3 *n* _{2} *⊗ J* _{2} 0

### ⎡

### ⎣

### ⎢ ⎢

### ⎢ ⎢

### ⎤

### ⎦

### ⎥ ⎥

### ⎥ ⎥

### ∈! ^{(n} ^{1} ^{n} ^{2} ^{)×(n} ^{1} ^{n} ^{2} ^{)}

^{(n}

^{n}

^{n}

### For FCC lattice

**Power method**

### Let be the eigenpairs of A where is linearly independent

### For any nonzero vector u

### Since , we have

### If for i >1 and , then

### Given shift value

### 9

### ( λ _{i} ^{, x} _{i} ^{) for i} *= 1,…,n*

_{i}

^{, x}

_{i}

^{) for i}

*x* _{1} , *…, x* _{n}

_{n}

*u* = α _{1} ^{x} _{1} +!+ α _{n} ^{x} _{n}

^{x}

_{n}

^{x}

_{n}

*A* ^{k} *x* _{i} = λ _{i} ^{k} ^{x} _{i}

^{k}

_{i}

_{i}

^{k}

^{x}

_{i}

*A* ^{k} *u* = α _{1} λ _{1} ^{k} ^{x} _{1} +!+ α _{n} λ _{n} ^{k} ^{x} _{n}

^{k}

^{k}

^{x}

_{n}

_{n}

^{k}

^{x}

_{n}

### λ _{1} > λ _{i} α _{1} ≠ 0

_{i}

### 1 λ _{1} ^{k} ^{A}

^{k}

^{A}

*k* *u* = α _{1} ^{x} _{1} + ( λ _{2} λ _{1} ^{)}

^{x}

*k* α _{2} ^{x} _{2} +!+ α _{n} ^{(} λ _{n} λ _{1} ^{)}

^{x}

_{n}

_{n}

*k* *x* _{n} → α _{1} ^{x} _{1} ^{ as k} → ∞

_{n}

^{x}

^{ as k}

*(A* − σ ^{I )} ^{−1}

^{I )}

### { } ^{k} ^{u} ^{=} ^{α} ^{1} { ^{(} ^{λ} ^{1} ^{−} ^{σ} ^{)} ^{−1} } ^{k} ^{x} ^{1} ^{+!+} ^{α} ^{n} { ^{(} ^{λ} ^{n} ^{−} ^{σ} ^{)} ^{−1} } ^{k} ^{x} ^{n}

^{k}

^{u}

^{k}

^{x}

^{n}

^{n}

^{k}

^{x}

^{n}

**Solving**

### Use shift-and-invert Lanczos method

### In each iteration of shift-and-invert Lanczos method, we need to solve

### How to efficiently solve this linear system?

### 10

*A* − λ ^{B}

^{B}

### ( ) ^{x} ^{= 0}

^{x}

*(A* − σ ^{B)y} *= b*

^{B)y}

**Solving linear system**

### 11

*(A* − σ ^{B)y} *= b*

^{B)y}

**Solve**

### Direct method (Gaussian elimination)

### Iterative method

### Matrix vector multiplication with Preconditioner M

### sol = bicgstabl(coef_mtx, rhs, tol, maxit,

### @(x)SSOR_prec(x, diag_coef_mtx, lower_L));

### 12

*(A* − σ ^{B)y} *= b*

^{B)y}

*y* *= (A −* σ ^{B) \ b}

^{B) \ b}

*A* − σ ^{B}

^{B}

**Solve**

### Direct method (Gaussian elimination)

### Iterative method

### Matrix vector multiplication with Preconditioner M

### sol = bicgstabl(coef_mtx, rhs, tol, maxit,

### @(x)SSOR_prec(x, diag_coef_mtx, lower_L));

### 12

*(A* − σ ^{B)y} *= b*

^{B)y}

*y* *= (A −* σ ^{B) \ b}

^{B) \ b}

*A* − σ ^{B}

^{B}

### Demo performance

**Eigen-decomp. of C** _{1} **, C** _{2} **, C** _{3} ** for SC lattice**

_{1}

_{2}

_{3}

### Define

### Define unitary matrix T as

### Then it holds that

### 13

*D* _{a,m} *= diag 1,e* ( ^{θ} ^{a ,m} , *!,e* ^{(m} ^{−1)} ^{θ} ^{a ,m} ) ^{,} ^{Λ} ^{a,m} ^{= diag e} ( ^{θ} ^{m ,1} ^{+} ^{θ} ^{a ,m} ^{−1 ! e} ^{θ} ^{m ,m} ^{+} ^{θ} ^{a ,m} ^{−1} ) ^{,}

_{a,m}

^{a ,m}^{(m}

^{a ,m}

^{a,m}^{= diag e}

^{m ,1}

^{a ,m}^{−1 ! e}

^{m ,m}

^{a ,m}*U* _{m} =

_{m}

### 1 1 ! 1

*e* ^{θ} ^{m ,1} *e* ^{θ} ^{m ,2} ! 1

^{m ,1}

^{m ,2}

### ! ! !

*e* ^{(m−1)} ^{θ} ^{m ,1} *e* ^{(m−1)} ^{θ} ^{m ,2} ! 1

^{(m−1)}

^{m ,1}

^{(m−1)}

^{m ,2}

### ⎡

### ⎣

### ⎢ ⎢

### ⎢ ⎢

### ⎤

### ⎦

### ⎥ ⎥

### ⎥ ⎥

### ∈! ^{m×m} , θ _{a,m} = *ı2* π ^{k} **⋅a**

^{m×m}

_{a,m}^{k}

*m* , θ _{m,i} = *ı2* π ^{i} *m*

_{m,i}

^{i}

*T* = 1

*n* *D* _{a}

_{a}

### 3 *,n* _{3} *⊗ D* _{a} _{2} _{,n} _{2} *⊗ D* _{a} _{1} _{,n} _{1}

_{a}

_{,n}

_{a}

_{,n}

### ( ) ( ^{U} ^{n} ^{3} ^{⊗U} ^{n} ^{2} ^{⊗U} ^{n} ^{1} )

^{U}

^{n}

^{⊗U}

^{n}

^{⊗U}

^{n}

*C* _{1} *T* = δ _{x} ^{−1} ^{T I} ( _{n} _{3} *⊗ I* _{n} _{2} ⊗ Λ _{a} _{1} _{,n} _{1} ) ^{≡ T Λ} ^{1} ^{,}

_{x}

^{T I}

_{n}

_{n}

_{a}

_{,n}

^{≡ T Λ}

*C* _{2} *T* = δ _{y} ^{−1} ^{T I} ( _{n} _{3} ⊗ Λ _{a} _{2} _{,n} _{2} *⊗ I* _{n} _{1} ) ^{≡ T Λ} ^{2} ^{,}

_{y}

^{T I}

_{n}

_{a}

_{,n}

_{n}

^{≡ T Λ}

*C* _{3} *T* = δ _{z} ^{−1} ^{T} ( Λ _{a} _{3} _{,n} _{3} *⊗ I* _{n} _{2} *⊗ I* _{n} _{1} ) ^{≡ T Λ} ^{3}

_{z}

^{T}

_{a}

_{,n}

_{n}

_{n}

^{≡ T Λ}

employing projections onto the range space of the discrete
matrix. (However, it is noted that the problem of many zeros
does not occur in the PWE as in each plane wave component
the number of polarizations is chosen to be the independent
two rather three.) (iii) As the third issue of importance, it is
of great interest to develop a fast algorithm for the
eigenvalue problem. An inverse method, accelerated by
multigrid technique with use of projection is proposed for
this purpose. The method exploits the sparsity of the matrix
for the eigenvalue problem in the finite diﬀerence formu-
lation. Because of the above mentioned diﬃculties, the
present method is a nontrivial extension of a similar method,
recently developed by the authors^{24)}for computing photonic
band structures in two dimensions.

In this study, we compute the band structures for three
types of photonic structures. The first one is a modified
simple cubic lattice consisting of dielectric spheres on the
lattice sites, each connected to its six nearest neighbors by
thin circular cylinders, which was proposed by Biswas et
al.^{25)}Figure 1 shows the modified simple cubic lattice. It is
noted that the original simple cubic structure comprising a
lattice of rods has been fabricated recently with advanced
silicon processing techniques.^{26)}The second one is the
tetragonal square spiral structure comprising a lattice of
circular or square cylinders, which was proposed by Toader
and John,^{27,28)}as shown in Fig. 2. Spiral structure was
discussed previously by Chutinan and Noda.^{29)}The square
spiral structure is arranged to connect the lattice points of
diamond structure with specific order, and is amenable to the
current technique of fabrication GLAD (GLancing Angle
Deposition) as discussed in refs. 30 and 31. As a third
example, we propose a diamond structure that has sp^{3}-like
configuration, composed of dielectric spheres with connect-
ing spheroids, as shown in Fig. 3. Diamond structures are
known to have large band gaps between relatively lower
branches either in diamond network or inverse diamond
structure.^{12,28)}In the present study, the spheroids, instead of
circular cylinders, take the positions of ‘‘valence bonds’’ to
imitate the sp^{3}structure of the electrons of diamond atoms.

Recently, submicron diamond-lattice photonic crystals have been successfully produced by two-photon laser nanofabri-

cation (photopolymerization).^{32)}

The order of presentation of the paper is organized as
follows. In §2, we show how to correctly formulate the finite
diﬀerence method for the double curl operator of the
photonic eigenvalue problem. In §3, we develop the
numerical method (inverse iteration with the full multigrid
acceleration) and present the fast algorithm, in which two
alternative methods of projection are proposed to avoid the
necessity of deflating zeros). In §4, we first present
numerical results that illustrate the eﬃciency of the
presently developed method. Then, the band structures are
computed for the modified simple cubic lattice, the tetrag-
onal square spiral structure (direct and inverse structure) and
the diamond structure with sp^{3}-like configuration. Finally,
concluding remarks with a summary of results are drawn in

§5.

2. Basic Equations and Finite Diﬀerence Formulation The electromagnetic waves propagating in the photonic crystals are well described by Maxwell’s equations. For linear isotropic and frequency-independent dielectric mate- rials with permeability close to one, the time-harmonic Fig. 1. Modified simple cubic lattice comprising dielectric spheres and

connecting thin circular cylinders.^{25)} Fig. 2. Tetragonal square spiral structure comprising circular cylin-
ders.^{27,28)}

Fig. 3. diamond structure with sp^{3}-like configuration comprising dielec-
tric spheres and connecting spheroids.

728 J. Phys. Soc. Jpn., Vol. 73, No. 3, March, 2004 R. L. CHERNet al.

**Eigen-decomp. of C** _{1} **, C** _{2} **, C** _{3} ** for FCC lattice**

_{1}

_{2}

_{3}

### Define

### Define unitary matrix T as

### Then it holds that

### 14

**x** _{i} *= D* _{x} *U* _{n} _{1} **(:,i), y** _{i, j} *= D* _{y,i} *U* _{n} _{2} *(:, j)*

_{i}

_{x}

_{n}

**(:,i), y**

_{i, j}

_{y,i}_{n}

*T* = 1

*n* ⎡ *T* _{1} *T* _{2} *! T* _{n} _{1}

_{n}

### ⎣⎢ ⎤

### ⎦⎥ ∈! ^{n} ^{×n} *, T* _{i} *= T* ⎡ _{i,1} *T* _{i,2} *! T* _{i,n} _{2}

^{n}

^{×n}

_{i}

_{i,1}

_{i,2}

_{i,n}

### ⎣⎢ ⎤

### ⎦⎥ ∈! ^{n} ^{×(n} ^{2} ^{n} ^{3} ^{)} , *T* _{i, j} *= D* _{z,i} _{+ j} *U* _{n}

^{n}

^{×(n}

^{n}

_{i, j}

_{z,i}_{+ j}

_{n}

### ( 3 ) ^{⊗ y} ^{(} ^{i, j} ^{⊗ x} ^{i} ^{)}

^{⊗ y}

^{i, j}

^{⊗ x}

^{i}

*C* _{1} *T* *= T Λ* ( _{n} _{1} *⊗ I* _{n} _{2} _{n} _{3} ) ^{≡ T Λ} ^{1} ^{,}

_{n}

_{n}

_{n}

^{≡ T Λ}

*C* _{2} *T* *= T ⊕* ( ( _{i=1} ^{n} ^{1} Λ _{i,n} _{2} ) ^{⊗ I} ^{n} ^{3} ) ^{≡ T Λ} ^{2} ^{,}

_{i=1}

^{n}

_{i,n}

^{⊗ I}

^{n}

^{≡ T Λ}

*C* _{3} *T* *= T ⊕* ( _{i=1} ^{n} ^{1} ⊕ ^{n} _{j} _{=1} ^{2} Λ _{i, j,n} _{3} ) ^{≡ T Λ} ^{3}

_{i=1}

^{n}

^{n}

_{j}

_{i, j,n}

^{≡ T Λ}

### ψ _{x} = *ı2* π ^{k} **⋅a** _{1}

_{x}

^{k}

*n* _{1} , *D* _{x} *= diag 1,e* ( ^{ψ}

_{x}

^{x}### , *!,e* ^{(n}

^{(n}

^{1}

^{−1)} ^{ψ}

^{x}### ) ^{,}

### ψ _{y,i} = *ı2* π

_{y,i}*n* _{2} **k** **⋅ a** _{2} − **a** _{1} 2

### ⎛ ⎝⎜ ⎞

### ⎠⎟ − *i* 2

### ⎧ ⎨

### ⎩

### ⎫ ⎬

### ⎭ , *D* _{y,i} *= diag 1,e* ( ^{ψ}

_{y,i}

^{y ,i}### , *!,e* ^{(n}

^{(n}

^{2}

^{−1)} ^{ψ}

^{y ,i}### ) ^{,}

### ψ _{z,i+ j} = *ı2* π

_{z,i+ j}*n* _{3} **k** **⋅ a** _{3} − **a** _{1} **+ a** _{2} 3

### ⎛ ⎝⎜ ⎞

### ⎠⎟ −

*i* *+ j* 3

### ⎧ ⎨

### ⎩

### ⎫ ⎬

### ⎭ *, D* _{z,i+ j} *= diag 1,e* ( ^{ψ}

_{z,i+ j}

^{y ,i+ j}### , *!,e* ^{(n}

^{(n}

^{3}

^{−1)} ^{ψ}

^{y ,i+ j}### )

^{24)}for computing photonic
band structures in two dimensions.

^{25)}Figure 1 shows the modified simple cubic lattice. It is
noted that the original simple cubic structure comprising a
lattice of rods has been fabricated recently with advanced
silicon processing techniques.^{26)} The second one is the
tetragonal square spiral structure comprising a lattice of
circular or square cylinders, which was proposed by Toader
and John,^{27,28)}as shown in Fig. 2. Spiral structure was
discussed previously by Chutinan and Noda.^{29)}The square
spiral structure is arranged to connect the lattice points of
diamond structure with specific order, and is amenable to the
current technique of fabrication GLAD (GLancing Angle
Deposition) as discussed in refs. 30 and 31. As a third
example, we propose a diamond structure that has sp^{3}-like
configuration, composed of dielectric spheres with connect-
ing spheroids, as shown in Fig. 3. Diamond structures are
known to have large band gaps between relatively lower
branches either in diamond network or inverse diamond
structure.^{12,28)}In the present study, the spheroids, instead of
circular cylinders, take the positions of ‘‘valence bonds’’ to
imitate the sp^{3}structure of the electrons of diamond atoms.

cation (photopolymerization).^{32)}

^{3}-like configuration. Finally,
concluding remarks with a summary of results are drawn in

§5.

connecting thin circular cylinders.^{25)} Fig. 2. Tetragonal square spiral structure comprising circular cylin-
ders.^{27,28)}

^{3}-like configuration comprising dielec-
tric spheres and connecting spheroids.

728 J. Phys. Soc. Jpn., Vol. 73, No. 3, March, 2004 R. L. CHERNet al.

**Eigen-decomp. of C** _{1} **, C** _{2} **, C** _{3} ** for FCC lattice**

_{1}

_{2}

_{3}

### Define

### Define unitary matrix T as

### Then it holds that

### 14

**x** _{i} *= D* _{x} *U* _{n} _{1} **(:,i), y** _{i, j} *= D* _{y,i} *U* _{n} _{2} *(:, j)*

_{i}

_{x}

_{n}

**(:,i), y**

_{i, j}

_{y,i}_{n}

*T* = 1

*n* ⎡ *T* _{1} *T* _{2} *! T* _{n} _{1}

_{n}

### ⎣⎢ ⎤

### ⎦⎥ ∈! ^{n} ^{×n} *, T* _{i} *= T* ⎡ _{i,1} *T* _{i,2} *! T* _{i,n} _{2}

^{n}

^{×n}

_{i}

_{i,1}

_{i,2}

_{i,n}

### ⎣⎢ ⎤

### ⎦⎥ ∈! ^{n} ^{×(n} ^{2} ^{n} ^{3} ^{)} , *T* _{i, j} *= D* _{z,i} _{+ j} *U* _{n}

^{n}

^{×(n}

^{n}

_{i, j}

_{z,i}_{+ j}

_{n}

### ( 3 ) ^{⊗ y} ^{(} ^{i, j} ^{⊗ x} ^{i} ^{)}

^{⊗ y}

^{i, j}

^{⊗ x}

^{i}

*C* _{1} *T* *= T Λ* ( _{n} _{1} *⊗ I* _{n} _{2} _{n} _{3} ) ^{≡ T Λ} ^{1} ^{,}

_{n}

_{n}

_{n}

^{≡ T Λ}

*C* _{2} *T* *= T ⊕* ( ( _{i=1} ^{n} ^{1} Λ _{i,n} _{2} ) ^{⊗ I} ^{n} ^{3} ) ^{≡ T Λ} ^{2} ^{,}

_{i=1}

^{n}

_{i,n}

^{⊗ I}

^{n}

^{≡ T Λ}

*C* _{3} *T* *= T ⊕* ( _{i=1} ^{n} ^{1} ⊕ ^{n} _{j} _{=1} ^{2} Λ _{i, j,n} _{3} ) ^{≡ T Λ} ^{3}

_{i=1}

^{n}

^{n}

_{j}

_{i, j,n}

^{≡ T Λ}

### ψ _{x} = *ı2* π ^{k} **⋅a** _{1}

_{x}

^{k}

*n* _{1} , *D* _{x} *= diag 1,e* ( ^{ψ}

_{x}

^{x}### , *!,e* ^{(n}

^{(n}

^{1}

^{−1)} ^{ψ}

^{x}### ) ^{,}

### ψ _{y,i} = *ı2* π

_{y,i}*n* _{2} **k** **⋅ a** _{2} − **a** _{1} 2

### ⎛ ⎝⎜ ⎞

### ⎠⎟ − *i* 2

### ⎧ ⎨

### ⎩

### ⎫ ⎬

### ⎭ , *D* _{y,i} *= diag 1,e* ( ^{ψ}

_{y,i}

^{y ,i}### , *!,e* ^{(n}

^{(n}

^{2}

^{−1)} ^{ψ}

^{y ,i}### ) ^{,}

### ψ _{z,i+ j} = *ı2* π

_{z,i+ j}*n* _{3} **k** **⋅ a** _{3} − **a** _{1} **+ a** _{2} 3

### ⎛ ⎝⎜ ⎞

### ⎠⎟ −

*i* *+ j* 3

### ⎧ ⎨

### ⎩

### ⎫ ⎬

### ⎭ *, D* _{z,i+ j} *= diag 1,e* ( ^{ψ}

_{z,i+ j}

^{y ,i+ j}### , *!,e* ^{(n}

^{(n}

^{3}

^{−1)} ^{ψ}

^{y ,i+ j}### )

^{24)}for computing photonic
band structures in two dimensions.

^{25)}Figure 1 shows the modified simple cubic lattice. It is
noted that the original simple cubic structure comprising a
lattice of rods has been fabricated recently with advanced
silicon processing techniques.^{26)} The second one is the
tetragonal square spiral structure comprising a lattice of
circular or square cylinders, which was proposed by Toader
and John,^{27,28)}as shown in Fig. 2. Spiral structure was
discussed previously by Chutinan and Noda.^{29)}The square
spiral structure is arranged to connect the lattice points of
diamond structure with specific order, and is amenable to the
current technique of fabrication GLAD (GLancing Angle
Deposition) as discussed in refs. 30 and 31. As a third
example, we propose a diamond structure that has sp^{3}-like
configuration, composed of dielectric spheres with connect-
ing spheroids, as shown in Fig. 3. Diamond structures are
known to have large band gaps between relatively lower
branches either in diamond network or inverse diamond
structure.^{12,28)}In the present study, the spheroids, instead of
circular cylinders, take the positions of ‘‘valence bonds’’ to
imitate the sp^{3}structure of the electrons of diamond atoms.

cation (photopolymerization).^{32)}

^{3}-like configuration. Finally,
concluding remarks with a summary of results are drawn in

§5.

connecting thin circular cylinders.^{25)} Fig. 2. Tetragonal square spiral structure comprising circular cylin-
ders.^{27,28)}

^{3}-like configuration comprising dielec-
tric spheres and connecting spheroids.

728 J. Phys. Soc. Jpn., Vol. 73, No. 3, March, 2004 R. L. CHERNet al.

### Demo performance

**CPU Times for T*p and Tq with FCC**

### 15

### Eigendecomposition of Double-Curl Operator 19

**0** **2** **4** **6** **8** **10**

**x 10** ^{7} **0**

^{7}

**1** **2** **3** **4** **5** **6** **7** **8** **9** **10**

** CPU times (sec.)**

** n** **Tq**

**T** ^{*} **p**

^{*}

### (a) n

**0** **1** **2** **3** **4** **5** **6** **7** **8**

**x 10** ^{8} **0**

^{8}

**1** **2** **3** **4** **5** **6** **7** **8** **9** **10**

** CPU times (sec.)**

** n log(n)** **Tq**

**T** ^{*} **p**

^{*}

### (b) n log(n) Fig. 6.1. CPU time for computing T ^{⇤} p and T q with various n.

### X U L G X W K

### 0 0.1 0.2 0.3 0.4 0.5 0.6 0.7 0.8 0.9

** frequency**

### Fig. 6.2. Band structure of the 3D photonic crystals with FCC lattice. The vectors k’s along the boundary of first Brillouin zone. The frequency ! = a p

### /(2⇡) is shown on the y-axis. The radius of the sphere is r = 0.12a and the connecting spheroid has minor axis length s = 0.11a.

### in MATLAB are used for the IPL and CG methods, respectively. The stopping tolerances for eigs and pcg are set to be 10 ^{4} ⇥ ✏/(2

### q 2

### x + _{y} ^{2} + _{z} ^{2} ) and ^{✏} 2

### x ⇥ 10 ^{3} , respectively, where ✏ ⇡ 2 ^{52} is the floating-point relative accuracy in MATLAB. The maximal number of Lanczos vectors for the restart in eigs is set to be 20. All computations are carried out on a workstation with two Intel Quad-Core Xeon X5687 3.6 GHz CPUs, 48 GB main memory, the RedHat Linux operation system, and IEEE double-precision floating-point arithmetic operations.

### Figure 6.1(a) shows the timing results for computing T ^{⇤} p and T q by Algorithm 2 and 3, respectively. The matrix size of T ranges from 884, 736 to 94, 818, 816. In particular, the dimension of T is ¯ n ^{3} _{j} , where ¯ n _{j} = 96 + 24j = n _{1} = n _{2} = n _{3} for j = 0, 1, . . . , 15. The average CPU time out of ten trials for each j is then plotted in the figure. We can see Algorithms 2 and 3 are extraordinarily efficient. They take less than 10 seconds to finish a T ^{⇤} p or T q matrix-vector multiplication even for the matrix T whose dimension is as large as 95 million. Figure 6.1(b) shows that the complexity of T ^{⇤} p and T q is O(n log(n)).

### Being equipped with these fast T ^{⇤} p or T q computational kernels, we evaluate how the IPL method (Algorithm 1) performs, in terms of CPU time and iteration numbers, to solve the eigenvalue problems for the band structure of the target photonic crystals.

**MATLAB** **T*p : ﬀt**

**Tq : iﬀt**

**Solving preconditioning linear system**

### 16

*(C* ^{∗} *C* − τ ^{I )y} ^{= d}

^{I )y}^{= d}

**Solving preconditioning linear system**

### 16

*C* ^{∗} *C* *= I* _{3} *⊗ G* ( ) ^{∗} *G* ^{− GG} ^{∗}

^{− GG}

*G* *= [C* _{1} ^{⊤} *,C* _{2} ^{⊤} *,C* _{3} ^{⊤} ] ^{⊤}

*C* =

### 0 *−C* _{3} *C* _{2} *C* _{3} 0 *−C* _{1}

*−C* _{2} *C* _{1} 0

### ⎡

### ⎣

### ⎢ ⎢

### ⎢

### ⎤

### ⎦

### ⎥ ⎥

### ⎥

*(C* ^{∗} *C* − τ ^{I )y} ^{= d}

^{I )y}^{= d}

**Solving preconditioning linear system**

### 16

*C* ^{∗} *C* *= I* _{3} *⊗ G* ( ) ^{∗} *G* ^{− GG} ^{∗}

^{− GG}

*G* *= [C* _{1} ^{⊤} *,C* _{2} ^{⊤} *,C* _{3} ^{⊤} ] ^{⊤}

*C* =

### 0 *−C* _{3} *C* _{2} *C* _{3} 0 *−C* _{1}

*−C* _{2} *C* _{1} 0

### ⎡

### ⎣

### ⎢ ⎢

### ⎢

### ⎤

### ⎦

### ⎥ ⎥

### ⎥

*(C* ^{∗} *C* − τ ^{I )y} ^{= d}

^{I )y}^{= d}

*I* _{3} *⊗ (G* ^{∗} *G)* − τ ^{I}

^{I}

### { } ^{y} ^{= d + GG} ^{∗} ^{y}

^{y}

^{= d + GG}^{y}

**Solving preconditioning linear system**

### 16

*C* ^{∗} *C* *= I* _{3} *⊗ G* ( ) ^{∗} *G* ^{− GG} ^{∗}

^{− GG}

*G* *= [C* _{1} ^{⊤} *,C* _{2} ^{⊤} *,C* _{3} ^{⊤} ] ^{⊤}

*C* =

### 0 *−C* _{3} *C* _{2} *C* _{3} 0 *−C* _{1}

*−C* _{2} *C* _{1} 0

### ⎡

### ⎣

### ⎢ ⎢

### ⎢

### ⎤

### ⎦

### ⎥ ⎥

### ⎥ *CG* = 0 _{GG} ^{∗} _{y} _{= −} _{τ} ^{−1} _{GG} ^{∗} _{d}

_{GG}

_{y}

_{GG}

_{d}

*(C* ^{∗} *C* − τ ^{I )y} ^{= d}

^{I )y}^{= d}

*I* _{3} *⊗ (G* ^{∗} *G)* − τ ^{I}

^{I}

### { } ^{y} ^{= d + GG} ^{∗} ^{y}

^{y}

^{= d + GG}^{y}

**Solving preconditioning linear system**

### 16

*C* ^{∗} *C* *= I* _{3} *⊗ G* ( ) ^{∗} *G* ^{− GG} ^{∗}

^{− GG}

*G* *= [C* _{1} ^{⊤} *,C* _{2} ^{⊤} *,C* _{3} ^{⊤} ] ^{⊤}

*C* =

### 0 *−C* _{3} *C* _{2} *C* _{3} 0 *−C* _{1}

*−C* _{2} *C* _{1} 0

### ⎡

### ⎣

### ⎢ ⎢

### ⎢

### ⎤

### ⎦

### ⎥ ⎥

### ⎥ *CG* = 0 _{GG} ^{∗} _{y} _{= −} _{τ} ^{−1} _{GG} ^{∗} _{d}

_{GG}

_{y}

_{GG}

_{d}

*(C* ^{∗} *C* − τ ^{I )y} ^{= d}

^{I )y}^{= d}

*I* _{3} *⊗ (G* ^{∗} *G)* − τ ^{I}

^{I}

### { } ^{y} ^{= d + GG} ^{∗} ^{y}

^{y}

^{= d + GG}^{y}

*I* _{3} *⊗ (G* ^{∗} *G)* − τ ^{I}

^{I}

### { } ^{y} ^{= d −} ^{τ} ^{−1} ^{GG} ^{∗} ^{d}

^{y}

^{= d −}

^{GG}

^{d}

**Solving preconditioning linear system**

### 16

*C* ^{∗} *C* *= I* _{3} *⊗ G* ( ) ^{∗} *G* ^{− GG} ^{∗}

^{− GG}

*G* *= [C* _{1} ^{⊤} *,C* _{2} ^{⊤} *,C* _{3} ^{⊤} ] ^{⊤}

*C* =

### 0 *−C* _{3} *C* _{2} *C* _{3} 0 *−C* _{1}

*−C* _{2} *C* _{1} 0

### ⎡

### ⎣

### ⎢ ⎢

### ⎢

### ⎤

### ⎦

### ⎥ ⎥

### ⎥ *CG* = 0 _{GG} ^{∗} _{y} _{= −} _{τ} ^{−1} _{GG} ^{∗} _{d}

_{GG}

_{y}

_{GG}

_{d}

*C* _{1} *T* *= T Λ* _{1} *, C* _{2} *T* *= T Λ* _{2} *, C* _{3} *T* *= T Λ* _{3} Λ _{q} = Λ _{1} ^{∗} Λ _{1} + Λ ^{∗} _{2} Λ _{2} + Λ ^{∗} _{3} Λ _{3}

_{q}

*(C* ^{∗} *C* − τ ^{I )y} ^{= d}

^{I )y}^{= d}

*I* _{3} *⊗ (G* ^{∗} *G)* − τ ^{I}

^{I}

### { } ^{y} ^{= d + GG} ^{∗} ^{y}

^{y}

^{= d + GG}^{y}

*I* _{3} *⊗ (G* ^{∗} *G)* − τ ^{I}

^{I}

### { } ^{y} ^{= d −} ^{τ} ^{−1} ^{GG} ^{∗} ^{d}

^{y}

^{= d −}

^{GG}

^{d}

**Solving preconditioning linear system**

### 16

*C* ^{∗} *C* *= I* _{3} *⊗ G* ( ) ^{∗} *G* ^{− GG} ^{∗}

^{− GG}

*G* *= [C* _{1} ^{⊤} *,C* _{2} ^{⊤} *,C* _{3} ^{⊤} ] ^{⊤}

*C* =

### 0 *−C* _{3} *C* _{2} *C* _{3} 0 *−C* _{1}

*−C* _{2} *C* _{1} 0

### ⎡

### ⎣

### ⎢ ⎢

### ⎢

### ⎤

### ⎦

### ⎥ ⎥

### ⎥ *CG* = 0 _{GG} ^{∗} _{y} _{= −} _{τ} ^{−1} _{GG} ^{∗} _{d}

_{GG}

_{y}

_{GG}

_{d}

*C* _{1} *T* *= T Λ* _{1} *, C* _{2} *T* *= T Λ* _{2} *, C* _{3} *T* *= T Λ* _{3} Λ _{q} = Λ _{1} ^{∗} Λ _{1} + Λ ^{∗} _{2} Λ _{2} + Λ ^{∗} _{3} Λ _{3}

_{q}

*(C* ^{∗} *C* − τ ^{I )y} ^{= d}

^{I )y}^{= d}

*I* _{3} *⊗ (G* ^{∗} *G)* − τ ^{I}

^{I}

### { } ^{y} ^{= d + GG} ^{∗} ^{y}

^{y}

^{= d + GG}^{y}

*I* _{3} *⊗ (G* ^{∗} *G)* − τ ^{I}

^{I}

### { } ^{y} ^{= d −} ^{τ} ^{−1} ^{GG} ^{∗} ^{d}

^{y}

^{= d −}

^{GG}

^{d}

*I* _{3} ⊗ Λ _{q} − τ ^{I}

_{q}

^{I}

### ( ) ^{y!} ^{= I −} ^{τ} ^{−1} ^{Λ} ^{Λ} ^{1} ^{2}

^{y!}

^{= I −}

### Λ _{3}

### ⎡

### ⎣

### ⎢ ⎢

### ⎢

### ⎤

### ⎦

### ⎥ ⎥

### ⎥ ⎡ Λ _{1} ^{*} Λ _{2} ^{*} Λ ^{*} _{3}

### ⎣⎢ ⎤

### ⎦⎥

### ⎛

### ⎝

### ⎜ ⎜

### ⎜

### ⎞

### ⎠

### ⎟ ⎟

### ⎟ ( *I* _{3} *⊗T* ) ^{*} ^{d, y} ^{= I} ( ^{3} ^{⊗T} ) ^{y!}

^{d, y}

^{= I}

^{⊗T}

^{y!}

**Solving preconditioning linear system**

### 16

*C* ^{∗} *C* *= I* _{3} *⊗ G* ( ) ^{∗} *G* ^{− GG} ^{∗}

^{− GG}

*G* *= [C* _{1} ^{⊤} *,C* _{2} ^{⊤} *,C* _{3} ^{⊤} ] ^{⊤}

*C* =

### 0 *−C* _{3} *C* _{2} *C* _{3} 0 *−C* _{1}

*−C* _{2} *C* _{1} 0

### ⎡

### ⎣

### ⎢ ⎢

### ⎢

### ⎤

### ⎦

### ⎥ ⎥

### ⎥ *CG* = 0 _{GG} ^{∗} _{y} _{= −} _{τ} ^{−1} _{GG} ^{∗} _{d}

_{GG}

_{y}

_{GG}

_{d}

*C* _{1} *T* *= T Λ* _{1} *, C* _{2} *T* *= T Λ* _{2} *, C* _{3} *T* *= T Λ* _{3} Λ _{q} = Λ _{1} ^{∗} Λ _{1} + Λ ^{∗} _{2} Λ _{2} + Λ ^{∗} _{3} Λ _{3}

_{q}

*(C* ^{∗} *C* − τ ^{I )y} ^{= d}

^{I )y}^{= d}

*I* _{3} *⊗ (G* ^{∗} *G)* − τ ^{I}

^{I}

### { } ^{y} ^{= d + GG} ^{∗} ^{y}

^{y}

^{= d + GG}^{y}

*I* _{3} *⊗ (G* ^{∗} *G)* − τ ^{I}

^{I}

### { } ^{y} ^{= d −} ^{τ} ^{−1} ^{GG} ^{∗} ^{d}

^{y}

^{= d −}

^{GG}

^{d}

*I* _{3} ⊗ Λ _{q} − τ ^{I}

_{q}

^{I}

### ( ) ^{y!} ^{= I −} ^{τ} ^{−1} ^{Λ} ^{Λ} ^{1} ^{2}

^{y!}

^{= I −}

### Λ _{3}

### ⎡

### ⎣

### ⎢ ⎢

### ⎢

### ⎤

### ⎦

### ⎥ ⎥

### ⎥ ⎡ Λ _{1} ^{*} Λ _{2} ^{*} Λ ^{*} _{3}

### ⎣⎢ ⎤

### ⎦⎥

### ⎛

### ⎝

### ⎜ ⎜

### ⎜

### ⎞

### ⎠

### ⎟ ⎟

### ⎟ ( *I* _{3} *⊗T* ) ^{*} ^{d, y} ^{= I} ( ^{3} ^{⊗T} ) ^{y!}

^{d, y}

^{= I}

^{⊗T}

^{y!}

### Demo performance

**Preconditioner**

### Iterative solver with preconditioner M:

### sol = bicgstabl(coef_mtx, rhs, tol, maxit, @(vec)FFT_based_precond(vec, Lambda, tau, EigDecompDoubCurl_cell, fun_mtx_TH_prod_vec,

### fun_mtx_T_prod_vec));

### 17

*M* *= C* ^{∗} *C* − τ ^{I}

^{I}

**Preconditioner**

### Iterative solver with preconditioner M:

### sol = bicgstabl(coef_mtx, rhs, tol, maxit, @(vec)FFT_based_precond(vec, Lambda, tau, EigDecompDoubCurl_cell, fun_mtx_TH_prod_vec,

### fun_mtx_T_prod_vec));

### 17

*M* *= C* ^{∗} *C* − τ ^{I}

^{I}

### Since

### we have

### No need to compute the matrix-vector multiplication involving A:

### sol = bicgstabl(@(vec)mtx_prod_vec_shift_invert_LS(vec, tau, Lambda_new, EigDecompDoubCurl_cell, mtx_B_sigma, fun_mtx_TH_prod_vec,

### fun_mtx_T_prod_vec), rhs, tol, maxit);

*M* ^{−1} *(A* − σ ^{B)} *= M* ^{−1} *(A- * τ ^{I} + τ ^{I} − σ ^{B)} *= I + M* ^{−1} ( τ ^{I} − σ ^{B)}

^{B)}

^{I}

^{I}

^{B)}

^{I}

^{B)}

*I* *+ M* ^{−1} ( τ *I* − σ *B)*

### { } ^{y} ^{= M} ^{−1} ^{b}

^{y}

^{= M}

^{b}

**Challenge in Solving Linear System**

### SC lattice (dim = 46875)

### FCC lattice

### 18

**Null-space free **

**eigenvalue problem **

### 19

**Huge zero eigenvalues**

### Eigen-decomposition

### where

### is unitary and

### 20 *Q* _{0} *Q*

### ⎡ ⎣ ⎤

*⎦ := I* ( ^{3} *⊗T* ) ^{⎡} _{⎣} ^{Π} ^{0} ^{Π} ^{1} ^{⎤} _{⎦ ≡ I} ( ^{3} ^{⊗T} ) ^{Π} ^{Π} ^{0,2} ^{0,1} ^{Π} ^{Π} ^{1,3} ^{1,1} ^{Π} ^{Π} ^{1,4} ^{1,2}

_{⎦ ≡ I}

^{⊗T}

### Π _{0,3} Π _{1,5} Π _{1,6}

### ⎡

### ⎣

### ⎢ ⎢

### ⎢ ⎢

### ⎤

### ⎦

### ⎥ ⎥

### ⎥ ⎥

*Q* _{0} *Q*

### ⎡ ⎣ ⎤

### ⎦

### ∗ *A* ⎡ *Q* _{0} *Q*

### ⎣ ⎤

### ⎦ = diag 0,Λ ( ^{q} , Λ _{q} ) ^{≡ diag 0,Λ} ^{( )}

^{q}

_{q}

*Q* ^{∗} *AQ* = Λ

### Λ _{q} = Λ _{1} ^{∗} Λ _{1} + Λ ^{∗} _{2} Λ _{2} + Λ ^{∗} _{3} Λ _{3}

_{q}