| 111 | |
| 112 | template<typename MAT, typename SCAL> |
| 113 | void CalcSchurComplement(const MAT a, FlatMatrix<SCAL> s, |
| 114 | const BitArray& used, LocalHeap& lh) |
| 115 | { |
| 116 | // static Timer t("HCurlAMG CalcSchurComplement"); RegionTimer reg(t); |
| 117 | |
| 118 | if (s.Height() == 0) return; |
| 119 | if (s.Height() == a.Height()) |
| 120 | { |
| 121 | s = a; |
| 122 | return; |
| 123 | } |
| 124 | |
| 125 | HeapReset hr(lh); |
| 126 | |
| 127 | int n = a.Height(); |
| 128 | Array<int> used_dofs(n, lh); |
| 129 | Array<int> unused_dofs(n, lh); |
| 130 | used_dofs.SetSize(0); |
| 131 | unused_dofs.SetSize(0); |
| 132 | for (int i = 0; i < n; i++) |
| 133 | if (used[i]) |
| 134 | used_dofs.Append(i); |
| 135 | else |
| 136 | unused_dofs.Append(i); |
| 137 | |
| 138 | s = a.Rows(used_dofs).Cols(used_dofs); |
| 139 | FlatMatrix<SCAL> b1 = a.Rows(unused_dofs).Cols(used_dofs) | lh; |
| 140 | FlatMatrix<SCAL> b2 = a.Rows(used_dofs).Cols(unused_dofs) | lh; |
| 141 | FlatMatrix<SCAL> c = a.Rows(unused_dofs).Cols(unused_dofs) | lh; |
| 142 | FlatMatrix<SCAL> hb1 (b1.Height(), b1.Width(), lh); |
| 143 | |
| 144 | CalcInverse (c); |
| 145 | hb1 = c * b1; |
| 146 | s -= b2 * hb1; |
| 147 | } |
| 148 | |
| 149 | |
| 150 | |