FieldMatVec.hpp 10.7 KB
(std::abs(b)); }; return std::pow( Impl::accumulate(x, op), 1.0/p ); } Praetorius, Simon committed Nov 13, 2017 243 244 245 246 /** \ingroup vector_norms * \brief The infty-norm of a vector = max_i |x_i| = alias for \ref abs_max **/ template Praetorius, Simon committed Jan 02, 2018 247 auto infty_norm(FieldVector const& x) Praetorius, Simon committed Nov 13, 2017 248 249 250 251 { return abs_max(x); } Praetorius, Simon committed Jan 22, 2018 252 253 254 255 256 257 template auto infty_norm(FieldMatrix const& x) { return abs_max(x); } Praetorius, Simon committed Nov 13, 2017 258 259 260 261 // ---------------------------------------------------------------------------- /// The euklidean distance between two vectors = |lhs-rhs|_2 template Praetorius, Simon committed Jan 02, 2018 262 T distance(FieldVector const& lhs, FieldVector const& rhs) Praetorius, Simon committed Nov 13, 2017 263 264 265 { T result = 0; for (int i = 0; i < N; ++i) Praetorius, Simon committed Dec 20, 2017 266 result += Math::sqr(lhs[i] - rhs[i]); Praetorius, Simon committed Nov 13, 2017 267 268 269 270 271 272 273 return std::sqrt(result); } // ---------------------------------------------------------------------------- /// Outer product (vec1 * vec2^T) template Praetorius, Simon committed Jan 02, 2018 274 auto outer(FieldMatrix const& vec1, FieldMatrix const& vec2) Praetorius, Simon committed Nov 13, 2017 275 { Praetorius, Simon committed Jan 02, 2018 276 using result_type = FieldMatrix() * std::declval() ), N, M>; Praetorius, Simon committed Nov 13, 2017 277 278 279 280 281 282 283 284 285 286 result_type mat; for (int i = 0; i < N; ++i) for (int j = 0; j < M; ++j) mat[i][j] = vec1[i].dot(vec2[j]); return mat; } // ---------------------------------------------------------------------------- template Praetorius, Simon committed Jan 02, 2018 287 T det(FieldMatrix const& /*mat*/) Praetorius, Simon committed Nov 13, 2017 288 289 290 291 292 293 { return 0; } /// Determinant of a 1x1 matrix template Praetorius, Simon committed Jan 02, 2018 294 T det(FieldMatrix const& mat) Praetorius, Simon committed Nov 13, 2017 295 296 297 298 299 300 { return mat[0][0]; } /// Determinant of a 2x2 matrix template Praetorius, Simon committed Jan 02, 2018 301 T det(FieldMatrix const& mat) Praetorius, Simon committed Nov 13, 2017 302 303 304 305 306 307 { return mat[0][0]*mat[1][1] - mat[0][1]*mat[1][0]; } /// Determinant of a 3x3 matrix template Praetorius, Simon committed Jan 02, 2018 308 T det(FieldMatrix const& mat) Praetorius, Simon committed Nov 13, 2017 309 310 311 312 313 314 315 { return mat[0][0]*mat[1][1]*mat[2][2] + mat[0][1]*mat[1][2]*mat[2][0] + mat[0][2]*mat[1][0]*mat[2][1] - (mat[0][2]*mat[1][1]*mat[2][0] + mat[0][1]*mat[1][0]*mat[2][2] + mat[0][0]*mat[1][2]*mat[2][1]); } /// Determinant of a NxN matrix template Praetorius, Simon committed Jan 02, 2018 316 T det(FieldMatrix const& mat) Praetorius, Simon committed Nov 13, 2017 317 318 319 320 321 322 { return mat.determinant(); } /// Return the inverse of the matrix mat template Praetorius, Simon committed Jan 02, 2018 323 auto inv(FieldMatrix mat) Praetorius, Simon committed Nov 13, 2017 324 325 326 327 328 329 330 { mat.invert(); return mat; } /// Solve the linear system A*x = b template Praetorius, Simon committed Jan 02, 2018 331 void solve(FieldMatrix const& A, FieldVector& x, FieldVector const& b) Praetorius, Simon committed Nov 13, 2017 332 333 334 335 336 337 338 { A.solve(x, b); } /// Gramian determinant = sqrt( det( DT^T * DF ) ) template Praetorius, Simon committed Jan 02, 2018 339 T gramian(FieldMatrix const& DF) Praetorius, Simon committed Nov 13, 2017 340 341 342 343 344 345 346 { using std::sqrt; return sqrt( det(outer(DF, DF)) ); } /// Gramian determinant, specialization for 1 column matrices template Praetorius, Simon committed Jan 02, 2018 347 T gramian(FieldMatrix const& DF) Praetorius, Simon committed Nov 13, 2017 348 349 350 351 352 { using std::sqrt; return sqrt(dot(DF[0], DF[0])); } Praetorius, Simon committed Dec 20, 2017 353 template Praetorius, Simon committed Jan 02, 2018 354 FieldMatrix trans(FieldMatrix const& A) Praetorius, Simon committed Dec 20, 2017 355 { Praetorius, Simon committed Jan 02, 2018 356 FieldMatrix At; Praetorius, Simon committed Dec 20, 2017 357 358 359 360 361 362 363 364 365 for (int i = 0; i < M; ++i) for (int j = 0; j < N; ++j) At[j][i] = A[i][j]; return At; } template Praetorius, Simon committed Jan 02, 2018 366 FieldMatrix multiplies(FieldMatrix const& A, FieldMatrix const& B) Praetorius, Simon committed Dec 20, 2017 367 368 369 370 371 372 373 { return A.rightmultiplyany(B); } template Praetorius, Simon committed Jan 02, 2018 374 FieldMatrix multiplies_AtB(FieldMatrix const& A, FieldMatrix const& B) Praetorius, Simon committed Dec 20, 2017 375 { Praetorius, Simon committed Jan 02, 2018 376 FieldMatrix C; Praetorius, Simon committed Dec 20, 2017 377 378 379 380 381 382 383 384 385 386 387 388 for (int m = 0; m < M; ++m) { for (int n = 0; n < N; ++n) { C[m][n] = 0; for (int l = 0; l < L; ++l) C[m][n] += A[l][m] * B[n][l]; } } return C; } template Praetorius, Simon committed Jan 02, 2018 389 FieldMatrix multiplies_ABt(FieldMatrix const& A, FieldMatrix const& B) Praetorius, Simon committed Dec 20, 2017 390 { Praetorius, Simon committed Jan 02, 2018 391 FieldMatrix C; Praetorius, Simon committed Dec 20, 2017 392 393 394 395 396 397 398 399 400 401 402 403 for (int m = 0; m < M; ++m) { for (int n = 0; n < N; ++n) { C[m][n] = 0; for (int l = 0; l < L; ++l) C[m][n] += A[m][l] * B[n][l]; } } return C; } template Praetorius, Simon committed Jan 02, 2018 404 FieldMatrix& multiplies_ABt(FieldMatrix const& A, FieldMatrix const& B, FieldMatrix& C) Praetorius, Simon committed Dec 20, 2017 405 406 407 408 409 410 411 412 413 414 415 { for (int m = 0; m < M; ++m) { for (int n = 0; n < N; ++n) { C[m][n] = 0; for (int l = 0; l < L; ++l) C[m][n] += A[m][l] * B[n][l]; } } return C; } Praetorius, Simon committed Jan 20, 2018 416 template Praetorius, Simon committed Jan 22, 2018 417 FieldMatrix& multiplies_ABt(FieldMatrix const& A, Dune::DiagonalMatrix const& B, FieldMatrix& C) Praetorius, Simon committed Jan 20, 2018 418 419 420 421 422 423 424 425 426 { for (int m = 0; m < M; ++m) { for (int n = 0; n < N; ++n) { C[m][n] = A[m][n] * B.diagonal(n); } } return C; } Praetorius, Simon committed Nov 13, 2017 427 } // end namespace AMDiS