File: mat.html

package info (click to toggle)
petsc 3.23.1%2Bdfsg1-1exp1
links: PTS, VCS
area: main
in suites: experimental
size: 515,576 kB
sloc: ansic: 751,607; cpp: 51,542; python: 38,598; f90: 17,352; javascript: 3,493; makefile: 3,157; sh: 1,502; xml: 619; objc: 445; java: 13; csh: 1
file content (2069 lines) | stat: -rw-r--r-- 205,152 bytes

<!DOCTYPE html>


<html lang="en" data-content_root="../" >

  <head>
    <meta charset="utf-8" />
    <meta name="viewport" content="width=device-width, initial-scale=1.0" /><meta name="viewport" content="width=device-width, initial-scale=1" />

    <title>Matrices &#8212; PETSc 3.23.1 documentation</title>
  
  
  
  <script data-cfasync="false">
    document.documentElement.dataset.mode = localStorage.getItem("mode") || "";
    document.documentElement.dataset.theme = localStorage.getItem("theme") || "light";
  </script>
  
  <!-- Loaded before other Sphinx assets -->
  <link href="../_static/styles/theme.css?digest=bd9e20870c6007c4c509" rel="stylesheet" />
<link href="../_static/styles/bootstrap.css?digest=bd9e20870c6007c4c509" rel="stylesheet" />
<link href="../_static/styles/pydata-sphinx-theme.css?digest=bd9e20870c6007c4c509" rel="stylesheet" />

  
  <link href="../_static/vendor/fontawesome/6.5.1/css/all.min.css?digest=bd9e20870c6007c4c509" rel="stylesheet" />
  <link rel="preload" as="font" type="font/woff2" crossorigin href="../_static/vendor/fontawesome/6.5.1/webfonts/fa-solid-900.woff2" />
<link rel="preload" as="font" type="font/woff2" crossorigin href="../_static/vendor/fontawesome/6.5.1/webfonts/fa-brands-400.woff2" />
<link rel="preload" as="font" type="font/woff2" crossorigin href="../_static/vendor/fontawesome/6.5.1/webfonts/fa-regular-400.woff2" />

    <link rel="stylesheet" type="text/css" href="../_static/pygments.css?v=8f2a1f02" />
    <link rel="stylesheet" type="text/css" href="../_static/copybutton.css?v=76b2166b" />
    <link rel="stylesheet" type="text/css" href="../_static/sphinx-design.min.css?v=87e54e7c" />
    <link rel="stylesheet" type="text/css" href="https://cdn.jsdelivr.net/npm/katex@0.16.10/dist/katex.min.css" />
    <link rel="stylesheet" type="text/css" href="../_static/katex-math.css?v=91adb8b6" />
    <link rel="stylesheet" type="text/css" href="../_static/css/custom.css?v=dbe1606d" />
  
  <!-- Pre-loaded scripts that we'll load fully later -->
  <link rel="preload" as="script" href="../_static/scripts/bootstrap.js?digest=bd9e20870c6007c4c509" />
<link rel="preload" as="script" href="../_static/scripts/pydata-sphinx-theme.js?digest=bd9e20870c6007c4c509" />
  <script src="../_static/vendor/fontawesome/6.5.1/js/all.min.js?digest=bd9e20870c6007c4c509"></script>

    <script src="../_static/documentation_options.js?v=34da53a5"></script>
    <script src="../_static/doctools.js?v=9a2dae69"></script>
    <script src="../_static/sphinx_highlight.js?v=dc90522c"></script>
    <script src="../_static/clipboard.min.js?v=a7894cd8"></script>
    <script src="../_static/copybutton.js?v=a56c686a"></script>
    <script src="../_static/design-tabs.js?v=f930bc37"></script>
    <script src="../_static/katex.min.js?v=be8ff15f"></script>
    <script src="../_static/auto-render.min.js?v=ad136472"></script>
    <script src="../_static/katex_autorenderer.js?v=bebc588a"></script>
    <script>DOCUMENTATION_OPTIONS.pagename = 'manual/mat';</script>
    <link rel="icon" href="../_static/petsc_favicon.png"/>
    <link rel="index" title="Index" href="../genindex.html" />
    <link rel="search" title="Search" href="../search.html" />
    <link rel="next" title="KSP: Linear System Solvers" href="ksp.html" />
    <link rel="prev" title="Vectors and Parallel Data" href="vec.html" />
  <meta name="viewport" content="width=device-width, initial-scale=1"/>
  <meta name="docsearch:language" content="en"/>
    <meta name="docbuild:last-update" content="2025-04-30T13:10:40-0500 (v3.23.1)"/>
  </head>
  
  
  <body data-bs-spy="scroll" data-bs-target=".bd-toc-nav" data-offset="180" data-bs-root-margin="0px 0px -60%" data-default-mode="">

  
  
  <a id="pst-skip-link" class="skip-link" href="#main-content">Skip to main content</a>
  
  <div id="pst-scroll-pixel-helper"></div>

  
  <button type="button" class="btn rounded-pill" id="pst-back-to-top">
    <i class="fa-solid fa-arrow-up"></i>
    Back to top
  </button>

  
  <input type="checkbox"
          class="sidebar-toggle"
          name="__primary"
          id="__primary"/>
  <label class="overlay overlay-primary" for="__primary"></label>
  
  <input type="checkbox"
          class="sidebar-toggle"
          name="__secondary"
          id="__secondary"/>
  <label class="overlay overlay-secondary" for="__secondary"></label>
  
  <div class="search-button__wrapper">
    <div class="search-button__overlay"></div>
    <div class="search-button__search-container">
<form class="bd-search d-flex align-items-center"
      action="../search.html"
      method="get">
  <i class="fa-solid fa-magnifying-glass"></i>
  <input type="search"
         class="form-control"
         name="q"
         id="search-input"
         placeholder="Search the docs ..."
         aria-label="Search the docs ..."
         autocomplete="off"
         autocorrect="off"
         autocapitalize="off"
         spellcheck="false"/>
  <span class="search-button__kbd-shortcut"><kbd class="kbd-shortcut__modifier">Ctrl</kbd>+<kbd>K</kbd></span>
</form></div>
  </div>

  <header>
  
    <div class="bd-header navbar navbar-expand-lg bd-navbar">
<div class="bd-header__inner bd-page-width">
  <label class="sidebar-toggle primary-toggle" for="__primary">
    <span class="fa-solid fa-bars"></span>
  </label>
  
  
  <div class="col-lg-3 navbar-header-items__start">
    
      <div class="navbar-item">

  

<a class="navbar-brand logo" href="../index.html">
  
  
  
  
  
    
    
      
    
    
    <img src="../_static/PETSc-TAO_RGB.svg" class="logo__image only-light" alt="PETSc 3.23.1 documentation - Home"/>
    <script>document.write(`<img src="../_static/PETSc-TAO_RGB_white.svg" class="logo__image only-dark" alt="PETSc 3.23.1 documentation - Home"/>`);</script>
  
  
</a></div>
    
  </div>
  
  <div class="col-lg-9 navbar-header-items">
    
    <div class="me-auto navbar-header-items__center">
      
        <div class="navbar-item">
<nav class="navbar-nav">
  <ul class="bd-navbar-elements navbar-nav">
    
                    <li class="nav-item current active">
                      <a class="nav-link nav-internal" href="../overview/index.html">
                        Overview
                      </a>
                    </li>
                

                    <li class="nav-item">
                      <a class="nav-link nav-internal" href="../install/index.html">
                        Install
                      </a>
                    </li>
                

                    <li class="nav-item">
                      <a class="nav-link nav-internal" href="../tutorials/index.html">
                        Tutorials
                      </a>
                    </li>
                

                    <li class="nav-item">
                      <a class="nav-link nav-internal" href="index.html">
                        User-Guide
                      </a>
                    </li>
                

                    <li class="nav-item">
                      <a class="nav-link nav-internal" href="../manualpages/index.html">
                        C/Fortran API
                      </a>
                    </li>
                

                    <li class="nav-item">
                      <a class="nav-link nav-internal" href="../petsc4py/index.html">
                        petsc4py API
                      </a>
                    </li>
                

                    <li class="nav-item">
                      <a class="nav-link nav-internal" href="../faq/index.html">
                        FAQ
                      </a>
                    </li>
                

                    <li class="nav-item">
                      <a class="nav-link nav-internal" href="../community/index.html">
                        Community
                      </a>
                    </li>
                

                    <li class="nav-item">
                      <a class="nav-link nav-internal" href="../developers/index.html">
                        Developers
                      </a>
                    </li>
                

                    <li class="nav-item">
                      <a class="nav-link nav-internal" href="../miscellaneous/index.html">
                        Misc.
                      </a>
                    </li>
                
  </ul>
</nav></div>
      
    </div>
    
    
    <div class="navbar-header-items__end">
      
        <div class="navbar-item navbar-persistent--container">
          

 <script>
 document.write(`
   <button class="btn navbar-btn search-button-field search-button__button" title="Search" aria-label="Search" data-bs-placement="bottom" data-bs-toggle="tooltip">
    <i class="fa-solid fa-magnifying-glass"></i>
    <span class="search-button__default-text">Search</span>
    <span class="search-button__kbd-shortcut"><kbd class="kbd-shortcut__modifier">Ctrl</kbd>+<kbd class="kbd-shortcut__modifier">K</kbd></span>
   </button>
 `);
 </script>
        </div>
      
      
        <div class="navbar-item">

<script>
document.write(`
  <button class="btn btn-sm navbar-btn theme-switch-button" title="light/dark" aria-label="light/dark" data-bs-placement="bottom" data-bs-toggle="tooltip">
    <span class="theme-switch nav-link" data-mode="light"><i class="fa-solid fa-sun fa-lg"></i></span>
    <span class="theme-switch nav-link" data-mode="dark"><i class="fa-solid fa-moon fa-lg"></i></span>
    <span class="theme-switch nav-link" data-mode="auto"><i class="fa-solid fa-circle-half-stroke fa-lg"></i></span>
  </button>
`);
</script></div>
      
        <div class="navbar-item"><ul class="navbar-icon-links navbar-nav"
    aria-label="Icon Links">
        <li class="nav-item">
          
          
          
          
          
          
          
          
          <a href="https://gitlab.com/petsc/petsc" title="GitLab" class="nav-link" rel="noopener" target="_blank" data-bs-toggle="tooltip" data-bs-placement="bottom"><span><i class="fab fa-gitlab fa-lg" aria-hidden="true"></i></span>
            <span class="sr-only">GitLab</span></a>
        </li>
</ul></div>
      
    </div>
    
  </div>
  
  
    <div class="navbar-persistent--mobile">

 <script>
 document.write(`
   <button class="btn navbar-btn search-button-field search-button__button" title="Search" aria-label="Search" data-bs-placement="bottom" data-bs-toggle="tooltip">
    <i class="fa-solid fa-magnifying-glass"></i>
    <span class="search-button__default-text">Search</span>
    <span class="search-button__kbd-shortcut"><kbd class="kbd-shortcut__modifier">Ctrl</kbd>+<kbd class="kbd-shortcut__modifier">K</kbd></span>
   </button>
 `);
 </script>
    </div>
  

  
    <label class="sidebar-toggle secondary-toggle" for="__secondary" tabindex="0">
      <span class="fa-solid fa-outdent"></span>
    </label>
  
</div>

    </div>
  
  </header>

  <div class="bd-container">
    <div class="bd-container__inner bd-page-width">
      
      
      
      <div class="bd-sidebar-primary bd-sidebar">
        

  
  <div class="sidebar-header-items sidebar-primary__section">
    
    
      <div class="sidebar-header-items__center">
        
          <div class="navbar-item">
<nav class="navbar-nav">
  <ul class="bd-navbar-elements navbar-nav">
    
                    <li class="nav-item current active">
                      <a class="nav-link nav-internal" href="../overview/index.html">
                        Overview
                      </a>
                    </li>
                

                    <li class="nav-item">
                      <a class="nav-link nav-internal" href="../install/index.html">
                        Install
                      </a>
                    </li>
                

                    <li class="nav-item">
                      <a class="nav-link nav-internal" href="../tutorials/index.html">
                        Tutorials
                      </a>
                    </li>
                

                    <li class="nav-item">
                      <a class="nav-link nav-internal" href="index.html">
                        User-Guide
                      </a>
                    </li>
                

                    <li class="nav-item">
                      <a class="nav-link nav-internal" href="../manualpages/index.html">
                        C/Fortran API
                      </a>
                    </li>
                

                    <li class="nav-item">
                      <a class="nav-link nav-internal" href="../petsc4py/index.html">
                        petsc4py API
                      </a>
                    </li>
                

                    <li class="nav-item">
                      <a class="nav-link nav-internal" href="../faq/index.html">
                        FAQ
                      </a>
                    </li>
                

                    <li class="nav-item">
                      <a class="nav-link nav-internal" href="../community/index.html">
                        Community
                      </a>
                    </li>
                

                    <li class="nav-item">
                      <a class="nav-link nav-internal" href="../developers/index.html">
                        Developers
                      </a>
                    </li>
                

                    <li class="nav-item">
                      <a class="nav-link nav-internal" href="../miscellaneous/index.html">
                        Misc.
                      </a>
                    </li>
                
  </ul>
</nav></div>
        
      </div>
    
    
    
      <div class="sidebar-header-items__end">
        
          <div class="navbar-item">

<script>
document.write(`
  <button class="btn btn-sm navbar-btn theme-switch-button" title="light/dark" aria-label="light/dark" data-bs-placement="bottom" data-bs-toggle="tooltip">
    <span class="theme-switch nav-link" data-mode="light"><i class="fa-solid fa-sun fa-lg"></i></span>
    <span class="theme-switch nav-link" data-mode="dark"><i class="fa-solid fa-moon fa-lg"></i></span>
    <span class="theme-switch nav-link" data-mode="auto"><i class="fa-solid fa-circle-half-stroke fa-lg"></i></span>
  </button>
`);
</script></div>
        
          <div class="navbar-item"><ul class="navbar-icon-links navbar-nav"
    aria-label="Icon Links">
        <li class="nav-item">
          
          
          
          
          
          
          
          
          <a href="https://gitlab.com/petsc/petsc" title="GitLab" class="nav-link" rel="noopener" target="_blank" data-bs-toggle="tooltip" data-bs-placement="bottom"><span><i class="fab fa-gitlab fa-lg" aria-hidden="true"></i></span>
            <span class="sr-only">GitLab</span></a>
        </li>
</ul></div>
        
      </div>
    
  </div>
  
    <div class="sidebar-primary-items__start sidebar-primary__section">
        <div class="sidebar-primary-item">
<nav class="bd-docs-nav bd-links"
     aria-label="Section Navigation">
  <p class="bd-links__title" role="heading" aria-level="1">Section Navigation</p>
  <div class="bd-toc-item navbar-nav"><ul class="current nav bd-sidenav">
<li class="toctree-l1"><a class="reference internal" href="../overview/nutshell.html">PETSc in a nutshell</a></li>
<li class="toctree-l1"><a class="reference internal" href="../overview/features.html">Supported Systems</a></li>

<li class="toctree-l1"><a class="reference internal" href="../overview/gpu_roadmap.html">GPU Support Roadmap</a></li>
<li class="toctree-l1"><a class="reference internal" href="../overview/vector_table.html">Summary of Vector Types Available In PETSc</a></li>
<li class="toctree-l1"><a class="reference internal" href="../overview/matrix_table.html">Summary of Matrix Types Available In PETSc</a></li>
<li class="toctree-l1"><a class="reference internal" href="../overview/linear_solve_table.html">Summary of Sparse Linear Solvers Available In PETSc</a></li>
<li class="toctree-l1"><a class="reference internal" href="../overview/nonlinear_solve_table.html">Summary of Nonlinear Solvers Available In PETSc</a></li>
<li class="toctree-l1"><a class="reference internal" href="../overview/integrator_table.html">Summary of Time Integrators Available In PETSc</a></li>
<li class="toctree-l1"><a class="reference internal" href="../overview/tao_solve_table.html">Summary of Tao Solvers</a></li>
<li class="toctree-l1"><a class="reference internal" href="../overview/discrete_table.html">Summary of Discretization Management Systems</a></li>
<li class="toctree-l1"><a class="reference internal" href="../overview/plex_transform_table.html">Summary of Unstructured Mesh Transformations</a></li>
<li class="toctree-l1 current active has-children"><a class="reference internal" href="index.html">User-Guide</a><input checked="" class="toctree-checkbox" id="toctree-checkbox-1" name="toctree-checkbox-1" type="checkbox"/><label class="toctree-toggle" for="toctree-checkbox-1"><i class="fa-solid fa-chevron-down"></i></label><ul class="current">
<li class="toctree-l2 has-children"><a class="reference internal" href="introduction.html">Introduction to PETSc</a><input class="toctree-checkbox" id="toctree-checkbox-2" name="toctree-checkbox-2" type="checkbox"/><label class="toctree-toggle" for="toctree-checkbox-2"><i class="fa-solid fa-chevron-down"></i></label><ul>
<li class="toctree-l3"><a class="reference internal" href="about_this_manual.html">About This Manual</a></li>
<li class="toctree-l3"><a class="reference internal" href="getting_started.html">Getting Started</a></li>






</ul>
</li>
<li class="toctree-l2 current active has-children"><a class="reference internal" href="programming.html">The Solvers in PETSc/TAO</a><input checked="" class="toctree-checkbox" id="toctree-checkbox-3" name="toctree-checkbox-3" type="checkbox"/><label class="toctree-toggle" for="toctree-checkbox-3"><i class="fa-solid fa-chevron-down"></i></label><ul class="current">
<li class="toctree-l3"><a class="reference internal" href="vec.html">Vectors and Parallel Data</a></li>
<li class="toctree-l3 current active"><a class="current reference internal" href="#">Matrices</a></li>
<li class="toctree-l3"><a class="reference internal" href="ksp.html">KSP: Linear System Solvers</a></li>
<li class="toctree-l3"><a class="reference internal" href="snes.html">SNES: Nonlinear Solvers</a></li>
<li class="toctree-l3"><a class="reference internal" href="ts.html">TS: Scalable ODE and DAE Solvers</a></li>

<li class="toctree-l3"><a class="reference internal" href="tao.html">TAO: Optimization Solvers</a></li>
</ul>
</li>
<li class="toctree-l2 has-children"><a class="reference internal" href="dm.html">DM: Interfacing Between Solvers and Models/Discretizations</a><input class="toctree-checkbox" id="toctree-checkbox-4" name="toctree-checkbox-4" type="checkbox"/><label class="toctree-toggle" for="toctree-checkbox-4"><i class="fa-solid fa-chevron-down"></i></label><ul>
<li class="toctree-l3"><a class="reference internal" href="dmbase.html">DM Basics</a></li>
<li class="toctree-l3"><a class="reference internal" href="section.html">PetscSection: Connecting Grids to Data</a></li>
<li class="toctree-l3"><a class="reference internal" href="dmplex.html">DMPlex: Unstructured Grids</a></li>
<li class="toctree-l3"><a class="reference internal" href="dmstag.html">DMSTAG: Staggered, Structured Grid</a></li>
<li class="toctree-l3"><a class="reference internal" href="dmnetwork.html">Networks</a></li>
<li class="toctree-l3"><a class="reference internal" href="dt.html">PetscDT: Discretization Technology in PETSc</a></li>
<li class="toctree-l3"><a class="reference internal" href="fe.html">PetscFE: Finite Element Infrastructure in PETSc</a></li>
</ul>
</li>
<li class="toctree-l2 has-children"><a class="reference internal" href="additional.html">Additional Information</a><input class="toctree-checkbox" id="toctree-checkbox-5" name="toctree-checkbox-5" type="checkbox"/><label class="toctree-toggle" for="toctree-checkbox-5"><i class="fa-solid fa-chevron-down"></i></label><ul>
<li class="toctree-l3"><a class="reference internal" href="fortran.html">PETSc for Fortran Users</a></li>
<li class="toctree-l3"><a class="reference internal" href="versionchecking.html">Checking the PETSc version</a></li>
<li class="toctree-l3"><a class="reference internal" href="matlab.html">Using MATLAB with PETSc</a></li>
<li class="toctree-l3"><a class="reference internal" href="profiling.html">Profiling</a></li>
<li class="toctree-l3"><a class="reference internal" href="performance.html">Hints for Performance Tuning</a></li>
<li class="toctree-l3"><a class="reference internal" href="streams.html">STREAMS: Example Study</a></li>
<li class="toctree-l3"><a class="reference internal" href="blas-lapack.html">The Use of BLAS and LAPACK in PETSc and external libraries</a></li>
<li class="toctree-l3"><a class="reference internal" href="other.html">Other PETSc Features</a></li>

<li class="toctree-l3"><a class="reference internal" href="advanced.html">Advanced Features of Matrices and Solvers</a></li>
<li class="toctree-l3"><a class="reference internal" href="tests.html">Running PETSc Tests</a></li>
</ul>
</li>
</ul>
</li>
<li class="toctree-l1 has-children"><a class="reference internal" href="../manualpages/index.html">C/Fortran API</a><input class="toctree-checkbox" id="toctree-checkbox-6" name="toctree-checkbox-6" type="checkbox"/><label class="toctree-toggle" for="toctree-checkbox-6"><i class="fa-solid fa-chevron-down"></i></label><ul>
<li class="toctree-l2 has-children"><a class="reference internal" href="../manualpages/Vector.html">Vectors and Index Sets</a><input class="toctree-checkbox" id="toctree-checkbox-7" name="toctree-checkbox-7" type="checkbox"/><label class="toctree-toggle" for="toctree-checkbox-7"><i class="fa-solid fa-chevron-down"></i></label><ul>
<li class="toctree-l3"><a class="reference internal" href="../manualpages/Vec/index.html">Vector Operations (Vec)</a></li>
<li class="toctree-l3"><a class="reference internal" href="../manualpages/IS/index.html">Index sets (IS)</a></li>
</ul>
</li>
<li class="toctree-l2 has-children"><a class="reference internal" href="../manualpages/Matrix.html">Matrices and Matrix Operations</a><input class="toctree-checkbox" id="toctree-checkbox-8" name="toctree-checkbox-8" type="checkbox"/><label class="toctree-toggle" for="toctree-checkbox-8"><i class="fa-solid fa-chevron-down"></i></label><ul>
<li class="toctree-l3"><a class="reference internal" href="../manualpages/Mat/index.html">Matrix Operations (Mat)</a></li>
<li class="toctree-l3"><a class="reference internal" href="../manualpages/MatGraphOperations/index.html">Matrix colorings (MatColoring), orderings (MatOrdering), partitionings (MatPartitioning), and coarsening (MatCoarsen)</a></li>
<li class="toctree-l3"><a class="reference internal" href="../manualpages/MatFD/index.html">Finite difference computation of Jacobians (MatFD)</a></li>
</ul>
</li>
<li class="toctree-l2 has-children"><a class="reference internal" href="../manualpages/DataLayout.html">Data Layout and Communication</a><input class="toctree-checkbox" id="toctree-checkbox-9" name="toctree-checkbox-9" type="checkbox"/><label class="toctree-toggle" for="toctree-checkbox-9"><i class="fa-solid fa-chevron-down"></i></label><ul>
<li class="toctree-l3"><a class="reference internal" href="../manualpages/PetscSF/index.html">Star Forest Communication (PetscSF)</a></li>
<li class="toctree-l3"><a class="reference internal" href="../manualpages/PetscSection/index.html">Section Data Layout (PetscSection)</a></li>
<li class="toctree-l3"><a class="reference internal" href="../manualpages/AO/index.html">Application Orderings (AO)</a></li>
</ul>
</li>
<li class="toctree-l2 has-children"><a class="reference internal" href="../manualpages/DataManagement.html">Data Management between Vec and Mat, and Distributed Mesh Data Structures</a><input class="toctree-checkbox" id="toctree-checkbox-10" name="toctree-checkbox-10" type="checkbox"/><label class="toctree-toggle" for="toctree-checkbox-10"><i class="fa-solid fa-chevron-down"></i></label><ul>
<li class="toctree-l3"><a class="reference internal" href="../manualpages/DM/index.html">Data Management (DM)</a></li>
<li class="toctree-l3"><a class="reference internal" href="../manualpages/DMDA/index.html">Structured Grids (DMDA)</a></li>
<li class="toctree-l3"><a class="reference internal" href="../manualpages/DMStag/index.html">Staggered, Structured Grids (DMSTAG)</a></li>
<li class="toctree-l3"><a class="reference internal" href="../manualpages/DMPlex/index.html">Unstructured Grids and Cell Complexes (DMPLEX)</a></li>
<li class="toctree-l3"><a class="reference internal" href="../manualpages/DMNetwork/index.html">Graphs and Networks (DMNETWORK)</a></li>
<li class="toctree-l3"><a class="reference internal" href="../manualpages/DMForest/index.html">A Forest of Trees and Structured Adaptive Refinement (DMFOREST)</a></li>
<li class="toctree-l3"><a class="reference internal" href="../manualpages/DMPatch/index.html">Sequences of parallel mesh patches (DMPATCH)</a></li>
<li class="toctree-l3"><a class="reference internal" href="../manualpages/DMSwarm/index.html">Particle Discretizations (DMSWARM)</a></li>
<li class="toctree-l3"><a class="reference internal" href="../manualpages/DMMOAB/index.html">MOAB Mesh Representation (DMMOAB)</a></li>
<li class="toctree-l3"><a class="reference internal" href="../manualpages/DMLabel/index.html">Selecting Parts of Meshes (DMLabel)</a></li>
<li class="toctree-l3"><a class="reference internal" href="../manualpages/DMPRODUCT/index.html">Tensor products of meshes (DMRODUCT)</a></li>
<li class="toctree-l3"><a class="reference internal" href="../manualpages/DMComposite/index.html">DMComposite</a></li>
</ul>
</li>
<li class="toctree-l2 has-children"><a class="reference internal" href="../manualpages/Discretization.html">Discretization and Function Spaces</a><input class="toctree-checkbox" id="toctree-checkbox-11" name="toctree-checkbox-11" type="checkbox"/><label class="toctree-toggle" for="toctree-checkbox-11"><i class="fa-solid fa-chevron-down"></i></label><ul>
<li class="toctree-l3"><a class="reference internal" href="../manualpages/DT/index.html">Discretization Technology and Quadrature (DT)</a></li>
<li class="toctree-l3"><a class="reference internal" href="../manualpages/SPACE/index.html">Function Spaces (PetscSpace)</a></li>
<li class="toctree-l3"><a class="reference internal" href="../manualpages/DUALSPACE/index.html">Dual Spaces (PetscDualSpace)</a></li>
<li class="toctree-l3"><a class="reference internal" href="../manualpages/FE/index.html">Finite Elements (PetscFE)</a></li>
<li class="toctree-l3"><a class="reference internal" href="../manualpages/FV/index.html">Finite Volumes (PetscFV)</a></li>
<li class="toctree-l3"><a class="reference internal" href="../manualpages/PF/index.html">Defining your own mathematical functions (PF)</a></li>
<li class="toctree-l3"><a class="reference internal" href="../manualpages/LANDAU/index.html">Landau Collision Operator</a></li>
</ul>
</li>
<li class="toctree-l2 has-children"><a class="reference internal" href="../manualpages/LinearSolvers.html">Linear Solvers and Preconditioners</a><input class="toctree-checkbox" id="toctree-checkbox-12" name="toctree-checkbox-12" type="checkbox"/><label class="toctree-toggle" for="toctree-checkbox-12"><i class="fa-solid fa-chevron-down"></i></label><ul>
<li class="toctree-l3"><a class="reference internal" href="../manualpages/KSP/index.html">Linear Solvers and Krylov Methods (KSP)</a></li>
<li class="toctree-l3"><a class="reference internal" href="../manualpages/PC/index.html">Preconditioners (PC)</a></li>
</ul>
</li>
<li class="toctree-l2 has-children"><a class="reference internal" href="../manualpages/NonlinearSolvers.html">Nonlinear Solvers</a><input class="toctree-checkbox" id="toctree-checkbox-13" name="toctree-checkbox-13" type="checkbox"/><label class="toctree-toggle" for="toctree-checkbox-13"><i class="fa-solid fa-chevron-down"></i></label><ul>
<li class="toctree-l3"><a class="reference internal" href="../manualpages/SNES/index.html">Nonlinear Solvers (SNES)</a></li>
<li class="toctree-l3"><a class="reference internal" href="../manualpages/SNESFAS/index.html">Full Approximation Scheme (FAS) nonlinear multigrid</a></li>
</ul>
</li>
<li class="toctree-l2 has-children"><a class="reference internal" href="../manualpages/Timestepping.html">Forward and Adjoint Timestepping</a><input class="toctree-checkbox" id="toctree-checkbox-14" name="toctree-checkbox-14" type="checkbox"/><label class="toctree-toggle" for="toctree-checkbox-14"><i class="fa-solid fa-chevron-down"></i></label><ul>
<li class="toctree-l3"><a class="reference internal" href="../manualpages/TS/index.html">Time Stepping ODE and DAE Solvers (TS)</a></li>
<li class="toctree-l3"><a class="reference internal" href="../manualpages/Sensitivity/index.html">Sensitivity Analysis for ODE and DAE</a></li>
<li class="toctree-l3"><a class="reference internal" href="../manualpages/Characteristic/index.html">Semi-Lagrangian Solves using the Method of Characteristics</a></li>
</ul>
</li>
<li class="toctree-l2 has-children"><a class="reference internal" href="../manualpages/Optimization.html">Optimization</a><input class="toctree-checkbox" id="toctree-checkbox-15" name="toctree-checkbox-15" type="checkbox"/><label class="toctree-toggle" for="toctree-checkbox-15"><i class="fa-solid fa-chevron-down"></i></label><ul>
<li class="toctree-l3"><a class="reference internal" href="../manualpages/Tao/index.html">Optimization Solvers (Tao)</a></li>
<li class="toctree-l3"><a class="reference internal" href="../manualpages/TaoLineSearch/index.html">Optimization Line Search (TaoLineSearch)</a></li>
</ul>
</li>
<li class="toctree-l2 has-children"><a class="reference internal" href="../manualpages/Visualization.html">Graphics and Visualization</a><input class="toctree-checkbox" id="toctree-checkbox-16" name="toctree-checkbox-16" type="checkbox"/><label class="toctree-toggle" for="toctree-checkbox-16"><i class="fa-solid fa-chevron-down"></i></label><ul>
<li class="toctree-l3"><a class="reference internal" href="../manualpages/Draw/index.html">Graphics (Draw)</a></li>
<li class="toctree-l3"><a class="reference internal" href="../manualpages/Viewer/index.html">Viewing Objects (Viewer)</a></li>
</ul>
</li>
<li class="toctree-l2 has-children"><a class="reference internal" href="../manualpages/System.html">System Routines, Profiling, Data Structures</a><input class="toctree-checkbox" id="toctree-checkbox-17" name="toctree-checkbox-17" type="checkbox"/><label class="toctree-toggle" for="toctree-checkbox-17"><i class="fa-solid fa-chevron-down"></i></label><ul>
<li class="toctree-l3"><a class="reference internal" href="../manualpages/Sys/index.html">PETSc Options, IO, and System Utilities</a></li>
<li class="toctree-l3"><a class="reference internal" href="../manualpages/PetscH/index.html">Hash Tables</a></li>
<li class="toctree-l3"><a class="reference internal" href="../manualpages/Log/index.html">Profiling and Logging</a></li>
<li class="toctree-l3"><a class="reference internal" href="../manualpages/Device/index.html">Device</a></li>
<li class="toctree-l3"><a class="reference internal" href="../manualpages/Matlab/index.html">Matlab</a></li>
<li class="toctree-l3"><a class="reference internal" href="../manualpages/Bag/index.html">Bag</a></li>
<li class="toctree-l3"><a class="reference internal" href="../manualpages/BM/index.html">Benchmarks (BM)</a></li>
</ul>
</li>
<li class="toctree-l2"><a class="reference internal" href="../changes/index.html">Changes for each release</a></li>
<li class="toctree-l2"><a class="reference internal" href="../manualpages/singleindex.html">Single Index of all PETSc Manual Pages</a></li>
</ul>
</li>
<li class="toctree-l1"><a class="reference internal" href="../changes/index.html">Changes for each release</a></li>
<li class="toctree-l1"><a class="reference internal" href="../manualpages/singleindex.html">Single Index of all PETSc Manual Pages</a></li>
<li class="toctree-l1 has-children"><a class="reference internal" href="../overview/previous_release_docs.html">Documentation for previous PETSc releases</a><input class="toctree-checkbox" id="toctree-checkbox-18" name="toctree-checkbox-18" type="checkbox"/><label class="toctree-toggle" for="toctree-checkbox-18"><i class="fa-solid fa-chevron-down"></i></label><ul>
<li class="toctree-l2"><a class="reference external" href="https://web.cels.anl.gov/projects/petsc/vault/petsc-3.22/docs"> 3.22</a></li>
<li class="toctree-l2"><a class="reference external" href="https://web.cels.anl.gov/projects/petsc/vault/petsc-3.21/docs"> 3.21</a></li>
<li class="toctree-l2"><a class="reference external" href="https://web.cels.anl.gov/projects/petsc/vault/petsc-3.20/docs"> 3.20</a></li>
<li class="toctree-l2"><a class="reference external" href="https://web.cels.anl.gov/projects/petsc/vault/petsc-3.19/docs"> 3.19</a></li>
<li class="toctree-l2"><a class="reference external" href="https://web.cels.anl.gov/projects/petsc/vault/petsc-3.18/docs"> 3.18</a></li>
<li class="toctree-l2"><a class="reference external" href="https://web.cels.anl.gov/projects/petsc/vault/petsc-3.17/docs"> 3.17</a></li>
<li class="toctree-l2"><a class="reference external" href="https://web.cels.anl.gov/projects/petsc/vault/petsc-3.16/docs"> 3.16</a></li>
<li class="toctree-l2"><a class="reference external" href="https://web.cels.anl.gov/projects/petsc/vault/petsc-3.15/docs"> 3.15</a></li>
<li class="toctree-l2"><a class="reference external" href="https://web.cels.anl.gov/projects/petsc/vault/petsc-3.14/docs"> 3.14</a></li>
<li class="toctree-l2"><a class="reference external" href="https://web.cels.anl.gov/projects/petsc/vault/petsc-3.13/docs"> 3.13</a></li>
<li class="toctree-l2"><a class="reference external" href="https://web.cels.anl.gov/projects/petsc/vault/petsc-3.12/docs"> 3.12</a></li>
<li class="toctree-l2"><a class="reference external" href="https://web.cels.anl.gov/projects/petsc/vault/petsc-3.11/docs"> 3.11</a></li>
<li class="toctree-l2"><a class="reference external" href="https://web.cels.anl.gov/projects/petsc/vault/petsc-3.10/docs"> 3.10</a></li>
<li class="toctree-l2"><a class="reference external" href="https://web.cels.anl.gov/projects/petsc/vault/petsc-3.9/docs"> 3.9</a></li>
<li class="toctree-l2"><a class="reference external" href="https://web.cels.anl.gov/projects/petsc/vault/petsc-3.8/docs"> 3.8</a></li>
<li class="toctree-l2"><a class="reference external" href="https://web.cels.anl.gov/projects/petsc/vault/petsc-3.7/docs"> 3.7</a></li>
<li class="toctree-l2"><a class="reference external" href="https://web.cels.anl.gov/projects/petsc/vault/petsc-3.6/docs"> 3.6</a></li>
<li class="toctree-l2"><a class="reference external" href="https://web.cels.anl.gov/projects/petsc/vault/petsc-3.5/docs"> 3.5</a></li>
<li class="toctree-l2"><a class="reference external" href="https://web.cels.anl.gov/projects/petsc/vault/petsc-3.4/docs"> 3.4</a></li>
<li class="toctree-l2"><a class="reference external" href="https://web.cels.anl.gov/projects/petsc/vault/petsc-3.3/docs"> 3.3</a></li>
</ul>
</li>
</ul>
</div>
</nav></div>
    </div>
  
  
  <div class="sidebar-primary-items__end sidebar-primary__section">
  </div>
  
  <div id="rtd-footer-container"></div>


      </div>
      
      <main id="main-content" class="bd-main">
        
        
          <div class="bd-content">
            <div class="bd-article-container">
              
              <div class="bd-header-article">
<div class="header-article-items header-article__inner">
  
    <div class="header-article-items__start">
      
        <div class="header-article-item">





<nav aria-label="Breadcrumb">
  <ul class="bd-breadcrumbs">
    
    <li class="breadcrumb-item breadcrumb-home">
      <a href="../index.html" class="nav-link" aria-label="Home">
        <i class="fa-solid fa-home"></i>
      </a>
    </li>
    
    <li class="breadcrumb-item"><a href="../overview/index.html" class="nav-link">Overview</a></li>
    
    
    <li class="breadcrumb-item"><i class="fa-solid fa-ellipsis"></i></li>
    
    
    <li class="breadcrumb-item"><a href="programming.html" class="nav-link">The Solvers in PETSc/TAO</a></li>
    
    <li class="breadcrumb-item active" aria-current="page">Matrices</li>
  </ul>
</nav>
</div>
      
    </div>
  
  
</div>
</div>
              
              
              
                
<div id="searchbox"></div>
                <article class="bd-article">
                  
  <section class="tex2jax_ignore mathjax_ignore" id="matrices">
<span id="ch-matrices"></span><h1>Matrices<a class="headerlink" href="#matrices" title="Link to this heading">#</a></h1>
<p>PETSc provides a variety of matrix implementations because no single
matrix format is appropriate for all problems. Currently, we support
dense storage and compressed sparse row storage (both sequential and
parallel versions) for CPU and GPU based matrices, as well as several specialized formats. Additional
specialized formats can be easily added.</p>
<p>This chapter describes the basics of using PETSc matrices in general
(regardless of the particular format chosen) and discusses tips for
efficient use of the several simple uniprocess and parallel matrix
types. The use of PETSc matrices involves the following actions: create
a particular type of matrix, insert values into it, process the matrix,
use the matrix for various computations, and finally destroy the matrix.
The application code does not need to know or care about the particular
storage formats of the matrices.</p>
<section id="creating-matrices">
<span id="sec-matcreate"></span><h2>Creating matrices<a class="headerlink" href="#creating-matrices" title="Link to this heading">#</a></h2>
<p>As with vectors, PETSc has APIs that allow the user to specify the exact details of the matrix
creation process but also <code class="docutils notranslate"><span class="pre"><a href="../manualpages/DM/DM.html">DM</a></span></code> based creation routines that handle most of the details automatically
for specific families of applications. This is done with</p>
<div class="highlight-c notranslate"><div class="highlight"><pre><span></span><span class="n"><a href="../manualpages/DM/DMCreateMatrix.html">DMCreateMatrix</a></span><span class="p">(</span><span class="n"><a href="../manualpages/DM/DM.html">DM</a></span><span class="w"> </span><span class="n">dm</span><span class="p">,</span><span class="n"><a href="../manualpages/Mat/Mat.html">Mat</a></span><span class="w"> </span><span class="o">*</span><span class="n">A</span><span class="p">)</span>
</pre></div>
</div>
<p>The type of matrix created can be controlled with either</p>
<div class="highlight-c notranslate"><div class="highlight"><pre><span></span><span class="n"><a href="../manualpages/DM/DMSetMatType.html">DMSetMatType</a></span><span class="p">(</span><span class="n"><a href="../manualpages/DM/DM.html">DM</a></span><span class="w"> </span><span class="n">dm</span><span class="p">,</span><span class="n"><a href="../manualpages/Mat/MatType.html">MatType</a></span><span class="w"> </span><span class="o">&lt;</span><span class="n"><a href="../manualpages/Mat/MATAIJ.html">MATAIJ</a></span><span class="w"> </span><span class="n">or</span><span class="w"> </span><span class="n"><a href="../manualpages/Mat/MATBAIJ.html">MATBAIJ</a></span><span class="w"> </span><span class="n">or</span><span class="w"> </span><span class="n"><a href="../manualpages/Mat/MATAIJCUSPARSE.html">MATAIJCUSPARSE</a></span><span class="w"> </span><span class="n">etc</span><span class="o">&gt;</span><span class="p">)</span>
</pre></div>
</div>
<p>or with</p>
<div class="highlight-c notranslate"><div class="highlight"><pre><span></span><span class="n"><a href="../manualpages/DM/DMSetFromOptions.html">DMSetFromOptions</a></span><span class="p">(</span><span class="n"><a href="../manualpages/DM/DM.html">DM</a></span><span class="w"> </span><span class="n">dm</span><span class="p">)</span>
</pre></div>
</div>
<p>and the options database option <code class="docutils notranslate"><span class="pre">-dm_mat_type</span> <span class="pre">&lt;aij</span> <span class="pre">or</span> <span class="pre">baij</span> <span class="pre">or</span> <span class="pre">aijcusparse</span> <span class="pre">etc&gt;</span></code> Matrices can be created for CPU usage, for GPU usage and for usage on
both the CPUs and GPUs.</p>
<p>The creation of <code class="docutils notranslate"><span class="pre"><a href="../manualpages/DM/DM.html">DM</a></span></code> objects is discussed in <a class="reference internal" href="vec.html#sec-struct"><span class="std std-ref">DMDA - Creating vectors for structured grids</span></a>, <a class="reference internal" href="vec.html#sec-unstruct"><span class="std std-ref">DMPLEX - Creating vectors for unstructured grids</span></a>, <a class="reference internal" href="vec.html#sec-network"><span class="std std-ref">DMNETWORK - Creating vectors for networks</span></a>.</p>
</section>
<section id="low-level-matrix-creation-routines">
<h2>Low-level matrix creation routines<a class="headerlink" href="#low-level-matrix-creation-routines" title="Link to this heading">#</a></h2>
<p>When using a <code class="docutils notranslate"><span class="pre"><a href="../manualpages/DM/DM.html">DM</a></span></code> is not practical for a particular application one can create matrices directly
using</p>
<div class="highlight-c notranslate"><div class="highlight"><pre><span></span><span class="n"><a href="../manualpages/Mat/MatCreate.html">MatCreate</a></span><span class="p">(</span><span class="n"><a href="../manualpages/Sys/MPI_Comm.html">MPI_Comm</a></span><span class="w"> </span><span class="n">comm</span><span class="p">,</span><span class="n"><a href="../manualpages/Mat/Mat.html">Mat</a></span><span class="w"> </span><span class="o">*</span><span class="n">A</span><span class="p">)</span>
<span class="n"><a href="../manualpages/Mat/MatSetSizes.html">MatSetSizes</a></span><span class="p">(</span><span class="n"><a href="../manualpages/Mat/Mat.html">Mat</a></span><span class="w"> </span><span class="n">A</span><span class="p">,</span><span class="n"><a href="../manualpages/Sys/PetscInt.html">PetscInt</a></span><span class="w"> </span><span class="n">m</span><span class="p">,</span><span class="n"><a href="../manualpages/Sys/PetscInt.html">PetscInt</a></span><span class="w"> </span><span class="n">n</span><span class="p">,</span><span class="n"><a href="../manualpages/Sys/PetscInt.html">PetscInt</a></span><span class="w"> </span><span class="n">M</span><span class="p">,</span><span class="n"><a href="../manualpages/Sys/PetscInt.html">PetscInt</a></span><span class="w"> </span><span class="n">N</span><span class="p">)</span>
</pre></div>
</div>
<p>This routine generates a sequential matrix when running one process and
a parallel matrix for two or more processes; the particular matrix
format is set by the user via options database commands. The user
specifies either the global matrix dimensions, given by <code class="docutils notranslate"><span class="pre">M</span></code> and <code class="docutils notranslate"><span class="pre">N</span></code>
or the local dimensions, given by <code class="docutils notranslate"><span class="pre">m</span></code> and <code class="docutils notranslate"><span class="pre">n</span></code> while PETSc completely
controls memory allocation. This routine facilitates switching among
various matrix types, for example, to determine the format that is most
efficient for a certain application. By default, <code class="docutils notranslate"><span class="pre"><a href="../manualpages/Mat/MatCreate.html">MatCreate</a>()</span></code> employs
the sparse AIJ format, which is discussed in detail in
<a class="reference internal" href="#sec-matsparse"><span class="std std-ref">Sparse Matrices</span></a>. See the manual pages for further
information about available matrix formats.</p>
</section>
<section id="assembling-putting-values-into-matrices">
<h2>Assembling (putting values into) matrices<a class="headerlink" href="#assembling-putting-values-into-matrices" title="Link to this heading">#</a></h2>
<p>To insert or add entries to a matrix on CPUs, one can call a variant of
<code class="docutils notranslate"><span class="pre"><a href="../manualpages/Mat/MatSetValues.html">MatSetValues</a>()</span></code>, either</p>
<div class="highlight-c notranslate"><div class="highlight"><pre><span></span><span class="n"><a href="../manualpages/Mat/MatSetValues.html">MatSetValues</a></span><span class="p">(</span><span class="n"><a href="../manualpages/Mat/Mat.html">Mat</a></span><span class="w"> </span><span class="n">A</span><span class="p">,</span><span class="n"><a href="../manualpages/Sys/PetscInt.html">PetscInt</a></span><span class="w"> </span><span class="n">m</span><span class="p">,</span><span class="k">const</span><span class="w"> </span><span class="n"><a href="../manualpages/Sys/PetscInt.html">PetscInt</a></span><span class="w"> </span><span class="n">idxm</span><span class="p">[],</span><span class="n"><a href="../manualpages/Sys/PetscInt.html">PetscInt</a></span><span class="w"> </span><span class="n">n</span><span class="p">,</span><span class="k">const</span><span class="w"> </span><span class="n"><a href="../manualpages/Sys/PetscInt.html">PetscInt</a></span><span class="w"> </span><span class="n">idxn</span><span class="p">[],</span><span class="k">const</span><span class="w"> </span><span class="n"><a href="../manualpages/Sys/PetscScalar.html">PetscScalar</a></span><span class="w"> </span><span class="n">values</span><span class="p">[],</span><span class="n"><a href="../manualpages/Sys/INSERT_VALUES.html">INSERT_VALUES</a></span><span class="p">);</span>
</pre></div>
</div>
<p>or</p>
<div class="highlight-c notranslate"><div class="highlight"><pre><span></span><span class="n"><a href="../manualpages/Mat/MatSetValues.html">MatSetValues</a></span><span class="p">(</span><span class="n"><a href="../manualpages/Mat/Mat.html">Mat</a></span><span class="w"> </span><span class="n">A</span><span class="p">,</span><span class="n"><a href="../manualpages/Sys/PetscInt.html">PetscInt</a></span><span class="w"> </span><span class="n">m</span><span class="p">,</span><span class="k">const</span><span class="w"> </span><span class="n"><a href="../manualpages/Sys/PetscInt.html">PetscInt</a></span><span class="w"> </span><span class="n">idxm</span><span class="p">[],</span><span class="n"><a href="../manualpages/Sys/PetscInt.html">PetscInt</a></span><span class="w"> </span><span class="n">n</span><span class="p">,</span><span class="k">const</span><span class="w"> </span><span class="n"><a href="../manualpages/Sys/PetscInt.html">PetscInt</a></span><span class="w"> </span><span class="n">idxn</span><span class="p">[],</span><span class="k">const</span><span class="w"> </span><span class="n"><a href="../manualpages/Sys/PetscScalar.html">PetscScalar</a></span><span class="w"> </span><span class="n">values</span><span class="p">[],</span><span class="n"><a href="../manualpages/Sys/ADD_VALUES.html">ADD_VALUES</a></span><span class="p">);</span>
</pre></div>
</div>
<p>This routine inserts or adds a logically dense subblock of dimension
<code class="docutils notranslate"><span class="pre">m*n</span></code> into the matrix. The integer indices <code class="docutils notranslate"><span class="pre">idxm</span></code> and <code class="docutils notranslate"><span class="pre">idxn</span></code>,
respectively, indicate the global row and column numbers to be inserted.
<code class="docutils notranslate"><span class="pre"><a href="../manualpages/Mat/MatSetValues.html">MatSetValues</a>()</span></code> uses the standard C convention, where the row and
column matrix indices begin with zero <em>regardless of the programming language
employed</em>. The array <code class="docutils notranslate"><span class="pre">values</span></code> is logically two-dimensional, containing
the values that are to be inserted. By default the values are given in
row major order, which is the opposite of the Fortran convention,
meaning that the value to be put in row <code class="docutils notranslate"><span class="pre">idxm[i]</span></code> and column
<code class="docutils notranslate"><span class="pre">idxn[j]</span></code> is located in <code class="docutils notranslate"><span class="pre">values[i*n+j]</span></code>. To allow the insertion of
values in column major order, one can call the command</p>
<div class="highlight-c notranslate"><div class="highlight"><pre><span></span><span class="n"><a href="../manualpages/Mat/MatSetOption.html">MatSetOption</a></span><span class="p">(</span><span class="n"><a href="../manualpages/Mat/Mat.html">Mat</a></span><span class="w"> </span><span class="n">A</span><span class="p">,</span><span class="n"><a href="../manualpages/Mat/MatOption.html">MAT_ROW_ORIENTED</a></span><span class="p">,</span><span class="n"><a href="../manualpages/Sys/PETSC_FALSE.html">PETSC_FALSE</a></span><span class="p">);</span>
</pre></div>
</div>
<p><em>Warning</em>: Several of the sparse implementations do <em>not</em> currently
support the column-oriented option.</p>
<p>This notation should not be a mystery to anyone. For example, to insert
one matrix into another when using MATLAB, one uses the command
<code class="docutils notranslate"><span class="pre">A(im,in)</span> <span class="pre">=</span> <span class="pre">B;</span></code> where <code class="docutils notranslate"><span class="pre">im</span></code> and <code class="docutils notranslate"><span class="pre">in</span></code> contain the indices for the
rows and columns. This action is identical to the calls above to
<code class="docutils notranslate"><span class="pre"><a href="../manualpages/Mat/MatSetValues.html">MatSetValues</a>()</span></code>.</p>
<p>When using the block compressed sparse row matrix format (<code class="docutils notranslate"><span class="pre"><a href="../manualpages/Mat/MATSEQBAIJ.html">MATSEQBAIJ</a></span></code>
or <code class="docutils notranslate"><span class="pre"><a href="../manualpages/Mat/MATMPIBAIJ.html">MATMPIBAIJ</a></span></code>), one can insert elements more efficiently using the
block variant, <code class="docutils notranslate"><span class="pre"><a href="../manualpages/Mat/MatSetValuesBlocked.html">MatSetValuesBlocked</a>()</span></code> or
<code class="docutils notranslate"><span class="pre"><a href="../manualpages/Mat/MatSetValuesBlockedLocal.html">MatSetValuesBlockedLocal</a>()</span></code>.</p>
<p>The function <code class="docutils notranslate"><span class="pre"><a href="../manualpages/Mat/MatSetOption.html">MatSetOption</a>()</span></code> accepts several other inputs; see the
manual page for details.</p>
<p>After the matrix elements have been inserted or added into the matrix,
they must be processed (also called “assembled”) before they can be
used. The routines for matrix processing are</p>
<div class="highlight-c notranslate"><div class="highlight"><pre><span></span><span class="n"><a href="../manualpages/Mat/MatAssemblyBegin.html">MatAssemblyBegin</a></span><span class="p">(</span><span class="n"><a href="../manualpages/Mat/Mat.html">Mat</a></span><span class="w"> </span><span class="n">A</span><span class="p">,</span><span class="n"><a href="../manualpages/Mat/MatAssemblyType.html">MAT_FINAL_ASSEMBLY</a></span><span class="p">);</span>
<span class="n"><a href="../manualpages/Mat/MatAssemblyEnd.html">MatAssemblyEnd</a></span><span class="p">(</span><span class="n"><a href="../manualpages/Mat/Mat.html">Mat</a></span><span class="w"> </span><span class="n">A</span><span class="p">,</span><span class="n"><a href="../manualpages/Mat/MatAssemblyType.html">MAT_FINAL_ASSEMBLY</a></span><span class="p">);</span>
</pre></div>
</div>
<p>By placing other code between these two calls, the user can perform
computations while messages are in transit. Calls to <code class="docutils notranslate"><span class="pre"><a href="../manualpages/Mat/MatSetValues.html">MatSetValues</a>()</span></code>
with the <code class="docutils notranslate"><span class="pre"><a href="../manualpages/Sys/INSERT_VALUES.html">INSERT_VALUES</a></span></code> and <code class="docutils notranslate"><span class="pre"><a href="../manualpages/Sys/ADD_VALUES.html">ADD_VALUES</a></span></code> options <em>cannot</em> be mixed
without intervening calls to the assembly routines. For such
intermediate assembly calls the second routine argument typically should
be <code class="docutils notranslate"><span class="pre"><a href="../manualpages/Mat/MatAssemblyType.html">MAT_FLUSH_ASSEMBLY</a></span></code>, which omits some of the work of the full
assembly process. <code class="docutils notranslate"><span class="pre"><a href="../manualpages/Mat/MatAssemblyType.html">MAT_FINAL_ASSEMBLY</a></span></code> is required only in the last
matrix assembly before a matrix is used.</p>
<p>Even though one may insert values into PETSc matrices without regard to
which process eventually stores them, for efficiency reasons we usually
recommend generating most entries on the process where they are destined
to be stored. To help the application programmer with this task for
matrices that are distributed across the processes by ranges, the
routine</p>
<div class="highlight-c notranslate"><div class="highlight"><pre><span></span><span class="n"><a href="../manualpages/Mat/MatGetOwnershipRange.html">MatGetOwnershipRange</a></span><span class="p">(</span><span class="n"><a href="../manualpages/Mat/Mat.html">Mat</a></span><span class="w"> </span><span class="n">A</span><span class="p">,</span><span class="n"><a href="../manualpages/Sys/PetscInt.html">PetscInt</a></span><span class="w"> </span><span class="o">*</span><span class="n">first_row</span><span class="p">,</span><span class="n"><a href="../manualpages/Sys/PetscInt.html">PetscInt</a></span><span class="w"> </span><span class="o">*</span><span class="n">last_row</span><span class="p">);</span>
</pre></div>
</div>
<p>informs the user that all rows from <code class="docutils notranslate"><span class="pre">first_row</span></code> to <code class="docutils notranslate"><span class="pre">last_row-1</span></code>
(since the value returned in <code class="docutils notranslate"><span class="pre">last_row</span></code> is one more than the global
index of the last local row) will be stored on the local process.</p>
<p>If the <code class="docutils notranslate"><span class="pre"><a href="../manualpages/Mat/Mat.html">Mat</a></span></code> was obtained from a <code class="docutils notranslate"><span class="pre"><a href="../manualpages/DM/DM.html">DM</a></span></code> with <code class="docutils notranslate"><span class="pre"><a href="../manualpages/DM/DMCreateMatrix.html">DMCreateMatrix</a>()</span></code>, then the range values are determined by the specific <code class="docutils notranslate"><span class="pre"><a href="../manualpages/DM/DM.html">DM</a></span></code>.
If the <code class="docutils notranslate"><span class="pre"><a href="../manualpages/Mat/Mat.html">Mat</a></span></code> was created directly, the range values are determined by the local sizes passed to <code class="docutils notranslate"><span class="pre"><a href="../manualpages/Mat/MatSetSizes.html">MatSetSizes</a>()</span></code> or <code class="docutils notranslate"><span class="pre"><a href="../manualpages/Mat/MatCreateAIJ.html">MatCreateAIJ</a>()</span></code> (and such low-level functions for other <code class="docutils notranslate"><span class="pre"><a href="../manualpages/Mat/MatType.html">MatType</a></span></code>).
If <code class="docutils notranslate"><span class="pre"><a href="../manualpages/Sys/PETSC_DECIDE.html">PETSC_DECIDE</a></span></code> was passed as the local size, then the vector uses default values for the range using <code class="docutils notranslate"><span class="pre"><a href="../manualpages/Sys/PetscSplitOwnership.html">PetscSplitOwnership</a>()</span></code>.
For certain <code class="docutils notranslate"><span class="pre"><a href="../manualpages/DM/DM.html">DM</a></span></code>, such as <code class="docutils notranslate"><span class="pre"><a href="../manualpages/DMDA/DMDA.html">DMDA</a></span></code>, it is better to use <code class="docutils notranslate"><span class="pre"><a href="../manualpages/DM/DM.html">DM</a></span></code> specific routines, such as <code class="docutils notranslate"><span class="pre"><a href="../manualpages/DMDA/DMDAGetGhostCorners.html">DMDAGetGhostCorners</a>()</span></code>, to determine
the local values in the matrix. See <a class="reference internal" href="#sec-matlayout"><span class="std std-ref">Matrix and Vector Layouts and Storage Locations</span></a> for full details on row and column layouts.</p>
<p>In the sparse matrix implementations, once the assembly routines have
been called, the matrices are compressed and can be used for
matrix-vector multiplication, etc. Any space for preallocated nonzeros
that was not filled by a call to <code class="docutils notranslate"><span class="pre"><a href="../manualpages/Mat/MatSetValues.html">MatSetValues</a>()</span></code> or a related routine
is compressed out by assembling with <code class="docutils notranslate"><span class="pre"><a href="../manualpages/Mat/MatAssemblyType.html">MAT_FINAL_ASSEMBLY</a></span></code>. If you
intend to use that extra space later, be sure to insert explicit zeros
before assembling with <code class="docutils notranslate"><span class="pre"><a href="../manualpages/Mat/MatAssemblyType.html">MAT_FINAL_ASSEMBLY</a></span></code> so the space will not be
compressed out. Once the matrix has been assembled, inserting new values
will be expensive since it will require copies and possible memory
allocation.</p>
<p>One may repeatedly assemble matrices that retain the same
nonzero pattern (such as within a nonlinear or time-dependent problem).
Where possible, data structures and communication
information will be reused (instead of regenerated) during successive
steps, thereby increasing efficiency. See
<a href="../src/ksp/ksp/tutorials/ex5.c.html">KSP Tutorial ex5</a>
for a simple example of solving two linear systems that use the same
matrix data structure.</p>
<p>For matrices associated with <code class="docutils notranslate"><span class="pre"><a href="../manualpages/DMDA/DMDA.html">DMDA</a></span></code> there is a higher-level interface for providing
the numerical values based on the concept of stencils. See the manual page of <code class="docutils notranslate"><span class="pre"><a href="../manualpages/Mat/MatSetValuesStencil.html">MatSetValuesStencil</a>()</span></code> for usage.</p>
<p>For GPUs the routines <code class="docutils notranslate"><span class="pre"><a href="../manualpages/Mat/MatSetPreallocationCOO.html">MatSetPreallocationCOO</a>()</span></code> and <code class="docutils notranslate"><span class="pre"><a href="../manualpages/Mat/MatSetValuesCOO.html">MatSetValuesCOO</a>()</span></code> should be used for efficient matrix assembly
instead of <code class="docutils notranslate"><span class="pre"><a href="../manualpages/Mat/MatSetValues.html">MatSetValues</a>()</span></code>.</p>
<p>We now introduce the various families of PETSc matrices. <code class="docutils notranslate"><span class="pre"><a href="../manualpages/DM/DMCreateMatrix.html">DMCreateMatrix</a>()</span></code> manages
the preallocation process (introduced below) automatically so many users do not need to
worry about the details of the preallocation process.</p>
<section id="matrix-and-vector-layouts-and-storage-locations">
<span id="sec-matlayout"></span><h3>Matrix and Vector Layouts and Storage Locations<a class="headerlink" href="#matrix-and-vector-layouts-and-storage-locations" title="Link to this heading">#</a></h3>
<p>The layout of PETSc matrices across MPI ranks is defined by two things</p>
<ul class="simple">
<li><p>the layout of the two compatible vectors in the computation of the matrix-vector product y = A * x and</p></li>
<li><p>the memory where various parts of the matrix are stored across the MPI ranks.</p></li>
</ul>
<p>PETSc vectors always have a contiguous range of vector entries stored on each MPI rank. The first rank has entries from 0 to <code class="docutils notranslate"><span class="pre">rend1</span></code> - 1, the
next rank has entries from <code class="docutils notranslate"><span class="pre">rend1</span></code> to <code class="docutils notranslate"><span class="pre">rend2</span></code> - 1, etc. Thus the ownership range on each rank is from <code class="docutils notranslate"><span class="pre">rstart</span></code> to <code class="docutils notranslate"><span class="pre">rend</span></code>, these values can be
obtained with <code class="docutils notranslate"><span class="pre"><a href="../manualpages/Vec/VecGetOwnershipRange.html">VecGetOwnershipRange</a></span></code>(<code class="docutils notranslate"><span class="pre"><a href="../manualpages/Vec/Vec.html">Vec</a></span></code> x, <code class="docutils notranslate"><span class="pre"><a href="../manualpages/Sys/PetscInt.html">PetscInt</a></span></code> * <code class="docutils notranslate"><span class="pre">rstart</span></code>, <code class="docutils notranslate"><span class="pre"><a href="../manualpages/Sys/PetscInt.html">PetscInt</a></span></code> * <code class="docutils notranslate"><span class="pre">rend</span></code>). Each PETSc <code class="docutils notranslate"><span class="pre"><a href="../manualpages/Vec/Vec.html">Vec</a></span></code> has a <code class="docutils notranslate"><span class="pre"><a href="../manualpages/IS/PetscLayout.html">PetscLayout</a></span></code> object that contains this information.</p>
<p>All PETSc matrices have two <code class="docutils notranslate"><span class="pre"><a href="../manualpages/IS/PetscLayout.html">PetscLayout</a></span></code>s, they define the vector layouts for y and x in the product, y = A * x. Their ownership range information
can be obtained with <code class="docutils notranslate"><span class="pre"><a href="../manualpages/Mat/MatGetOwnershipRange.html">MatGetOwnershipRange</a>()</span></code>, <code class="docutils notranslate"><span class="pre"><a href="../manualpages/Mat/MatGetOwnershipRangeColumn.html">MatGetOwnershipRangeColumn</a>()</span></code>, <code class="docutils notranslate"><span class="pre"><a href="../manualpages/Mat/MatGetOwnershipRanges.html">MatGetOwnershipRanges</a>()</span></code>, and <code class="docutils notranslate"><span class="pre"><a href="../manualpages/Mat/MatGetOwnershipRangesColumn.html">MatGetOwnershipRangesColumn</a>()</span></code>.
Note that <code class="docutils notranslate"><span class="pre"><a href="../manualpages/Mat/MatCreateVecs.html">MatCreateVecs</a>()</span></code> provides two vectors that have compatible layouts for the associated vector.</p>
<p>For most PETSc matrices, excluding <code class="docutils notranslate"><span class="pre"><a href="../manualpages/Mat/MATELEMENTAL.html">MATELEMENTAL</a></span></code> and <code class="docutils notranslate"><span class="pre"><a href="../manualpages/Mat/MATSCALAPACK.html">MATSCALAPACK</a></span></code>, the row ownership range obtained with <code class="docutils notranslate"><span class="pre"><a href="../manualpages/Mat/MatGetOwnershipRange.html">MatGetOwnershipRange</a>()</span></code> also defines
where the matrix entries are stored; the matrix entries for rows <code class="docutils notranslate"><span class="pre">rstart</span></code> to <code class="docutils notranslate"><span class="pre">rend</span> <span class="pre">-</span> <span class="pre">1</span></code> are stored on the corresponding MPI rank. For other matrices
the rank where each matrix entry is stored is more complicated; information about the storage locations can be obtained with <code class="docutils notranslate"><span class="pre"><a href="../manualpages/Mat/MatGetOwnershipIS.html">MatGetOwnershipIS</a>()</span></code>.
Note that for
most PETSc matrices the values returned by <code class="docutils notranslate"><span class="pre"><a href="../manualpages/Mat/MatGetOwnershipIS.html">MatGetOwnershipIS</a>()</span></code> are the same as those returned by <code class="docutils notranslate"><span class="pre"><a href="../manualpages/Mat/MatGetOwnershipRange.html">MatGetOwnershipRange</a>()</span></code> and
<code class="docutils notranslate"><span class="pre"><a href="../manualpages/Mat/MatGetOwnershipRangeColumn.html">MatGetOwnershipRangeColumn</a>()</span></code>.</p>
<p>The PETSc object <code class="docutils notranslate"><span class="pre"><a href="../manualpages/IS/PetscLayout.html">PetscLayout</a></span></code> contains the ownership information that is provided by <code class="docutils notranslate"><span class="pre"><a href="../manualpages/Vec/VecGetOwnershipRange.html">VecGetOwnershipRange</a>()</span></code> and with <code class="docutils notranslate"><span class="pre"><a href="../manualpages/Mat/MatGetOwnershipRange.html">MatGetOwnershipRange</a>()</span></code>, <code class="docutils notranslate"><span class="pre"><a href="../manualpages/Mat/MatGetOwnershipRangeColumn.html">MatGetOwnershipRangeColumn</a>()</span></code>. Each vector has one layout, which can be obtained with <code class="docutils notranslate"><span class="pre"><a href="../manualpages/Vec/VecGetLayout.html">VecGetLayout</a>()</span></code> and <code class="docutils notranslate"><span class="pre"><a href="../manualpages/Mat/MatGetLayouts.html">MatGetLayouts</a>()</span></code>. Layouts support the routines <code class="docutils notranslate"><span class="pre"><a href="../manualpages/IS/PetscLayoutGetLocalSize.html">PetscLayoutGetLocalSize</a>()</span></code>, <code class="docutils notranslate"><span class="pre"><a href="../manualpages/IS/PetscLayoutGetSize.html">PetscLayoutGetSize</a>()</span></code>, <code class="docutils notranslate"><span class="pre"><a href="../manualpages/IS/PetscLayoutGetBlockSize.html">PetscLayoutGetBlockSize</a>()</span></code>, <code class="docutils notranslate"><span class="pre"><a href="../manualpages/IS/PetscLayoutGetRanges.html">PetscLayoutGetRanges</a>()</span></code>, <code class="docutils notranslate"><span class="pre"><a href="../manualpages/IS/PetscLayoutCompare.html">PetscLayoutCompare</a>()</span></code> as well as a variety of creation routines. These are used by the <code class="docutils notranslate"><span class="pre"><a href="../manualpages/Vec/Vec.html">Vec</a></span></code> and <code class="docutils notranslate"><span class="pre"><a href="../manualpages/Mat/Mat.html">Mat</a></span></code> and so are rarely needed directly. Finally <code class="docutils notranslate"><span class="pre"><a href="../manualpages/Sys/PetscSplitOwnership.html">PetscSplitOwnership</a>()</span></code> is a utility routine that does the same splitting of ownership ranges as <code class="docutils notranslate"><span class="pre"><a href="../manualpages/IS/PetscLayout.html">PetscLayout</a></span></code>.</p>
</section>
<section id="sparse-matrices">
<span id="sec-matsparse"></span><h3>Sparse Matrices<a class="headerlink" href="#sparse-matrices" title="Link to this heading">#</a></h3>
<p>The default matrix representation within PETSc is the general sparse AIJ
format (also called the compressed sparse
row format, CSR). This section discusses tips for <em>efficiently</em> using
this matrix format for large-scale applications. Additional formats
(such as block compressed row and block symmetric storage, which are
generally much more efficient for problems with multiple degrees of
freedom per node) are discussed below. Beginning users need not concern
themselves initially with such details and may wish to proceed directly
to <a class="reference internal" href="#sec-matoptions"><span class="std std-ref">Basic Matrix Operations</span></a>. However, when an application code
progresses to the point of tuning for efficiency and/or generating
timing results, it is <em>crucial</em> to read this information.</p>
<section id="sequential-aij-sparse-matrices">
<h4>Sequential AIJ Sparse Matrices<a class="headerlink" href="#sequential-aij-sparse-matrices" title="Link to this heading">#</a></h4>
<p>In the PETSc AIJ matrix formats, we store the nonzero elements by rows,
along with an array of corresponding column numbers and an array of
pointers to the beginning of each row. Note that the diagonal matrix
entries are stored with the rest of the nonzeros (not separately).</p>
<p>To create a sequential AIJ sparse matrix, <code class="docutils notranslate"><span class="pre">A</span></code>, with <code class="docutils notranslate"><span class="pre">m</span></code> rows and
<code class="docutils notranslate"><span class="pre">n</span></code> columns, one uses the command</p>
<div class="highlight-c notranslate"><div class="highlight"><pre><span></span><span class="n"><a href="../manualpages/Mat/MatCreateSeqAIJ.html">MatCreateSeqAIJ</a></span><span class="p">(</span><span class="n"><a href="../manualpages/Sys/PETSC_COMM_SELF.html">PETSC_COMM_SELF</a></span><span class="p">,</span><span class="n"><a href="../manualpages/Sys/PetscInt.html">PetscInt</a></span><span class="w"> </span><span class="n">m</span><span class="p">,</span><span class="n"><a href="../manualpages/Sys/PetscInt.html">PetscInt</a></span><span class="w"> </span><span class="n">n</span><span class="p">,</span><span class="n"><a href="../manualpages/Sys/PetscInt.html">PetscInt</a></span><span class="w"> </span><span class="n">nz</span><span class="p">,</span><span class="n"><a href="../manualpages/Sys/PetscInt.html">PetscInt</a></span><span class="w"> </span><span class="o">*</span><span class="n">nnz</span><span class="p">,</span><span class="n"><a href="../manualpages/Mat/Mat.html">Mat</a></span><span class="w"> </span><span class="o">*</span><span class="n">A</span><span class="p">);</span>
</pre></div>
</div>
<p>where <code class="docutils notranslate"><span class="pre">nz</span></code> or <code class="docutils notranslate"><span class="pre">nnz</span></code> can be used to preallocate matrix memory, as
discussed below. The user can set <code class="docutils notranslate"><span class="pre">nz=0</span></code> and <code class="docutils notranslate"><span class="pre">nnz=NULL</span></code> for PETSc to
control all matrix memory allocation.</p>
<p>The sequential and parallel AIJ matrix storage formats by default employ
<em>i-nodes</em> (identical nodes) when possible. We search for consecutive
rows with the same nonzero structure, thereby reusing matrix information
for increased efficiency. Related options database keys are
<code class="docutils notranslate"><span class="pre">-mat_no_inode</span></code> (do not use i-nodes) and <code class="docutils notranslate"><span class="pre">-mat_inode_limit</span> <span class="pre">&lt;limit&gt;</span></code>
(set i-node limit (max limit=5)). Note that problems with a single degree
of freedom per grid node will automatically not use i-nodes.</p>
<p>The internal data representation for the AIJ formats employs zero-based
indexing.</p>
</section>
<section id="preallocation-of-memory-for-sequential-aij-sparse-matrices">
<h4>Preallocation of Memory for Sequential AIJ Sparse Matrices<a class="headerlink" href="#preallocation-of-memory-for-sequential-aij-sparse-matrices" title="Link to this heading">#</a></h4>
<p>The dynamic process of allocating new memory and copying from the old
storage to the new is <em>intrinsically very expensive</em>. Thus, to obtain
good performance when assembling an AIJ matrix, it is crucial to
preallocate the memory needed for the sparse matrix. The user has two
choices for preallocating matrix memory via <code class="docutils notranslate"><span class="pre"><a href="../manualpages/Mat/MatCreateSeqAIJ.html">MatCreateSeqAIJ</a>()</span></code>.</p>
<p>One can use the scalar <code class="docutils notranslate"><span class="pre">nz</span></code> to specify the expected number of nonzeros
for each row. This is generally fine if the number of nonzeros per row
is roughly the same throughout the matrix (or as a quick and easy first
step for preallocation). If one underestimates the actual number of
nonzeros in a given row, then during the assembly process PETSc will
automatically allocate additional needed space. However, this extra
memory allocation can slow the computation.</p>
<p>If different rows have very different numbers of nonzeros, one should
attempt to indicate (nearly) the exact number of elements intended for
the various rows with the optional array, <code class="docutils notranslate"><span class="pre">nnz</span></code> of length <code class="docutils notranslate"><span class="pre">m</span></code>, where
<code class="docutils notranslate"><span class="pre">m</span></code> is the number of rows, for example</p>
<div class="highlight-c notranslate"><div class="highlight"><pre><span></span><span class="n"><a href="../manualpages/Sys/PetscInt.html">PetscInt</a></span><span class="w"> </span><span class="n">nnz</span><span class="p">[</span><span class="n">m</span><span class="p">];</span>
<span class="n">nnz</span><span class="p">[</span><span class="mi">0</span><span class="p">]</span><span class="w"> </span><span class="o">=</span><span class="w"> </span><span class="o">&lt;</span><span class="n">nonzeros</span><span class="w"> </span><span class="n">in</span><span class="w"> </span><span class="n">row</span><span class="w"> </span><span class="mi">0</span><span class="o">&gt;</span>
<span class="n">nnz</span><span class="p">[</span><span class="mi">1</span><span class="p">]</span><span class="w"> </span><span class="o">=</span><span class="w"> </span><span class="o">&lt;</span><span class="n">nonzeros</span><span class="w"> </span><span class="n">in</span><span class="w"> </span><span class="n">row</span><span class="w"> </span><span class="mi">1</span><span class="o">&gt;</span>
<span class="p">....</span>
<span class="n">nnz</span><span class="p">[</span><span class="n">m</span><span class="mi">-1</span><span class="p">]</span><span class="w"> </span><span class="o">=</span><span class="w"> </span><span class="o">&lt;</span><span class="n">nonzeros</span><span class="w"> </span><span class="n">in</span><span class="w"> </span><span class="n">row</span><span class="w"> </span><span class="n">m</span><span class="mi">-1</span><span class="o">&gt;</span>
</pre></div>
</div>
<p>In this case, the assembly process will require no additional memory
allocations if the <code class="docutils notranslate"><span class="pre">nnz</span></code> estimates are correct. If, however, the
<code class="docutils notranslate"><span class="pre">nnz</span></code> estimates are incorrect, PETSc will automatically obtain the
additional needed space, at a slight loss of efficiency.</p>
<p>Using the array <code class="docutils notranslate"><span class="pre">nnz</span></code> to preallocate memory is especially important
for efficient matrix assembly if the number of nonzeros varies
considerably among the rows. One can generally set <code class="docutils notranslate"><span class="pre">nnz</span></code> either by
knowing in advance the problem structure (e.g., the stencil for finite
difference problems on a structured grid) or by precomputing the
information by using a segment of code similar to that for the regular
matrix assembly. The overhead of determining the <code class="docutils notranslate"><span class="pre">nnz</span></code> array will be
quite small compared with the overhead of the inherently expensive
<code class="docutils notranslate"><span class="pre">malloc</span></code>s and moves of data that are needed for dynamic allocation
during matrix assembly. Always guess high if an exact value is not known
(extra space is cheaper than too little).</p>
<p>Thus, when assembling a sparse matrix with very different numbers of
nonzeros in various rows, one could proceed as follows for finite
difference methods:</p>
<ol class="arabic simple">
<li><p>Allocate integer array <code class="docutils notranslate"><span class="pre">nnz</span></code>.</p></li>
<li><p>Loop over grid, counting the expected number of nonzeros for the
row(s) associated with the various grid points.</p></li>
<li><p>Create the sparse matrix via <code class="docutils notranslate"><span class="pre"><a href="../manualpages/Mat/MatCreateSeqAIJ.html">MatCreateSeqAIJ</a>()</span></code> or alternative.</p></li>
<li><p>Loop over the grid, generating matrix entries and inserting in matrix
via <code class="docutils notranslate"><span class="pre"><a href="../manualpages/Mat/MatSetValues.html">MatSetValues</a>()</span></code>.</p></li>
</ol>
<p>For (vertex-based) finite element type calculations, an analogous
procedure is as follows:</p>
<ol class="arabic simple">
<li><p>Allocate integer array <code class="docutils notranslate"><span class="pre">nnz</span></code>.</p></li>
<li><p>Loop over vertices, computing the number of neighbor vertices, which
determines the number of nonzeros for the corresponding matrix
row(s).</p></li>
<li><p>Create the sparse matrix via <code class="docutils notranslate"><span class="pre"><a href="../manualpages/Mat/MatCreateSeqAIJ.html">MatCreateSeqAIJ</a>()</span></code> or alternative.</p></li>
<li><p>Loop over elements, generating matrix entries and inserting in matrix
via <code class="docutils notranslate"><span class="pre"><a href="../manualpages/Mat/MatSetValues.html">MatSetValues</a>()</span></code>.</p></li>
</ol>
<p>The <code class="docutils notranslate"><span class="pre">-info</span></code> option causes the routines <code class="docutils notranslate"><span class="pre"><a href="../manualpages/Mat/MatAssemblyBegin.html">MatAssemblyBegin</a>()</span></code> and
<code class="docutils notranslate"><span class="pre"><a href="../manualpages/Mat/MatAssemblyEnd.html">MatAssemblyEnd</a>()</span></code> to print information about the success of the
preallocation. Consider the following example for the <code class="docutils notranslate"><span class="pre"><a href="../manualpages/Mat/MATSEQAIJ.html">MATSEQAIJ</a></span></code>
matrix format:</p>
<div class="highlight-c notranslate"><div class="highlight"><pre><span></span><span class="nl">MatAssemblyEnd_SeqAIJ</span><span class="p">:</span><span class="n">Matrix</span><span class="w"> </span><span class="n">size</span><span class="w"> </span><span class="mi">10</span><span class="w"> </span><span class="n">X</span><span class="w"> </span><span class="mi">10</span><span class="p">;</span><span class="w"> </span><span class="n">storage</span><span class="w"> </span><span class="n">space</span><span class="o">:</span><span class="mi">20</span><span class="w"> </span><span class="n">unneeded</span><span class="p">,</span><span class="w"> </span><span class="mi">100</span><span class="w"> </span><span class="n">used</span>
<span class="nl">MatAssemblyEnd_SeqAIJ</span><span class="p">:</span><span class="n">Number</span><span class="w"> </span><span class="n">of</span><span class="w"> </span><span class="n">mallocs</span><span class="w"> </span><span class="n">during</span><span class="w"> </span><span class="n"><a href="../manualpages/Mat/MatSetValues.html">MatSetValues</a></span><span class="w"> </span><span class="n">is</span><span class="w"> </span><span class="mi">0</span>
</pre></div>
</div>
<p>The first line indicates that the user preallocated 120 spaces but only
100 were used. The second line indicates that the user preallocated
enough space so that PETSc did not have to internally allocate
additional space (an expensive operation). In the next example the user
did not preallocate sufficient space, as indicated by the fact that the
number of mallocs is very large (bad for efficiency):</p>
<div class="highlight-c notranslate"><div class="highlight"><pre><span></span><span class="nl">MatAssemblyEnd_SeqAIJ</span><span class="p">:</span><span class="n">Matrix</span><span class="w"> </span><span class="n">size</span><span class="w"> </span><span class="mi">10</span><span class="w"> </span><span class="n">X</span><span class="w"> </span><span class="mi">10</span><span class="p">;</span><span class="w"> </span><span class="n">storage</span><span class="w"> </span><span class="n">space</span><span class="o">:</span><span class="mi">47</span><span class="w"> </span><span class="n">unneeded</span><span class="p">,</span><span class="w"> </span><span class="mi">1000</span><span class="w"> </span><span class="n">used</span>
<span class="nl">MatAssemblyEnd_SeqAIJ</span><span class="p">:</span><span class="n">Number</span><span class="w"> </span><span class="n">of</span><span class="w"> </span><span class="n">mallocs</span><span class="w"> </span><span class="n">during</span><span class="w"> </span><span class="n"><a href="../manualpages/Mat/MatSetValues.html">MatSetValues</a></span><span class="w"> </span><span class="n">is</span><span class="w"> </span><span class="mi">40000</span>
</pre></div>
</div>
<p>Although at first glance such procedures for determining the matrix
structure in advance may seem unusual, they are actually very efficient
because they alleviate the need for dynamic construction of the matrix
data structure, which can be very expensive.</p>
</section>
<section id="parallel-aij-sparse-matrices">
<h4>Parallel AIJ Sparse Matrices<a class="headerlink" href="#parallel-aij-sparse-matrices" title="Link to this heading">#</a></h4>
<p>Parallel sparse matrices with the AIJ format can be created with the
command</p>
<div class="highlight-c notranslate"><div class="highlight"><pre><span></span><span class="n"><a href="../manualpages/Mat/MatCreateAIJ.html">MatCreateAIJ</a></span><span class="p">(</span><span class="n"><a href="../manualpages/Sys/MPI_Comm.html">MPI_Comm</a></span><span class="w"> </span><span class="n">comm</span><span class="p">,</span><span class="n"><a href="../manualpages/Sys/PetscInt.html">PetscInt</a></span><span class="w"> </span><span class="n">m</span><span class="p">,</span><span class="n"><a href="../manualpages/Sys/PetscInt.html">PetscInt</a></span><span class="w"> </span><span class="n">n</span><span class="p">,</span><span class="n"><a href="../manualpages/Sys/PetscInt.html">PetscInt</a></span><span class="w"> </span><span class="n">M</span><span class="p">,</span><span class="n"><a href="../manualpages/Sys/PetscInt.html">PetscInt</a></span><span class="w"> </span><span class="n">N</span><span class="p">,</span><span class="n"><a href="../manualpages/Sys/PetscInt.html">PetscInt</a></span><span class="w"> </span><span class="n">d_nz</span><span class="p">,</span><span class="n"><a href="../manualpages/Sys/PetscInt.html">PetscInt</a></span><span class="w"> </span><span class="o">*</span><span class="n">d_nnz</span><span class="p">,</span><span class="w"> </span><span class="n"><a href="../manualpages/Sys/PetscInt.html">PetscInt</a></span><span class="w"> </span><span class="n">o_nz</span><span class="p">,</span><span class="n"><a href="../manualpages/Sys/PetscInt.html">PetscInt</a></span><span class="w"> </span><span class="o">*</span><span class="n">o_nnz</span><span class="p">,</span><span class="n"><a href="../manualpages/Mat/Mat.html">Mat</a></span><span class="w"> </span><span class="o">*</span><span class="n">A</span><span class="p">);</span>
</pre></div>
</div>
<p><code class="docutils notranslate"><span class="pre">A</span></code> is the newly created matrix, while the arguments <code class="docutils notranslate"><span class="pre">m</span></code>, <code class="docutils notranslate"><span class="pre">M</span></code>, and
<code class="docutils notranslate"><span class="pre">N</span></code>, indicate the number of local rows and the number of global rows
and columns, respectively. In the PETSc partitioning scheme, all the
matrix columns are local and <code class="docutils notranslate"><span class="pre">n</span></code> is the number of columns
corresponding to the local part of a parallel vector. Either the local or
global parameters can be replaced with <code class="docutils notranslate"><span class="pre"><a href="../manualpages/Sys/PETSC_DECIDE.html">PETSC_DECIDE</a></span></code>, so that PETSc
will determine them. The matrix is stored with a fixed number of rows on
each process, given by <code class="docutils notranslate"><span class="pre">m</span></code>, or determined by PETSc if <code class="docutils notranslate"><span class="pre">m</span></code> is
<code class="docutils notranslate"><span class="pre"><a href="../manualpages/Sys/PETSC_DECIDE.html">PETSC_DECIDE</a></span></code>.</p>
<p>If <code class="docutils notranslate"><span class="pre"><a href="../manualpages/Sys/PETSC_DECIDE.html">PETSC_DECIDE</a></span></code> is not used for the arguments <code class="docutils notranslate"><span class="pre">m</span></code> and <code class="docutils notranslate"><span class="pre">n</span></code>, then
the user must ensure that they are chosen to be compatible with the
vectors. To do this, one first considers the matrix-vector product
<span class="math">\(y = A x\)</span>. The <code class="docutils notranslate"><span class="pre">m</span></code> that is used in the matrix creation routine
<code class="docutils notranslate"><span class="pre"><a href="../manualpages/Mat/MatCreateAIJ.html">MatCreateAIJ</a>()</span></code> must match the local size used in the vector creation
routine <code class="docutils notranslate"><span class="pre"><a href="../manualpages/Vec/VecCreateMPI.html">VecCreateMPI</a>()</span></code> for <code class="docutils notranslate"><span class="pre">y</span></code>. Likewise, the <code class="docutils notranslate"><span class="pre">n</span></code> used must
match that used as the local size in <code class="docutils notranslate"><span class="pre"><a href="../manualpages/Vec/VecCreateMPI.html">VecCreateMPI</a>()</span></code> for <code class="docutils notranslate"><span class="pre">x</span></code>.</p>
<p>The user must set <code class="docutils notranslate"><span class="pre">d_nz=0</span></code>, <code class="docutils notranslate"><span class="pre">o_nz=0</span></code>, <code class="docutils notranslate"><span class="pre">d_nnz=</span></code>NULL, and
<code class="docutils notranslate"><span class="pre">o_nnz=NULL</span></code> for PETSc to control dynamic allocation of matrix memory
space. Analogous to <code class="docutils notranslate"><span class="pre">nz</span></code> and <code class="docutils notranslate"><span class="pre">nnz</span></code> for the routine
<code class="docutils notranslate"><span class="pre"><a href="../manualpages/Mat/MatCreateSeqAIJ.html">MatCreateSeqAIJ</a>()</span></code>, these arguments optionally specify nonzero
information for the diagonal (<code class="docutils notranslate"><span class="pre">d_nz</span></code> and <code class="docutils notranslate"><span class="pre">d_nnz</span></code>) and off-diagonal
(<code class="docutils notranslate"><span class="pre">o_nz</span></code> and <code class="docutils notranslate"><span class="pre">o_nnz</span></code>) parts of the matrix. For a square global
matrix, we define each process’s diagonal portion to be its local rows
and the corresponding columns (a square submatrix); each process’s
off-diagonal portion encompasses the remainder of the local matrix (a
rectangular submatrix). The rank in the MPI communicator determines the
absolute ordering of the blocks. That is, the process with rank 0 in the
communicator given to <code class="docutils notranslate"><span class="pre"><a href="../manualpages/Mat/MatCreateAIJ.html">MatCreateAIJ</a>()</span></code> contains the top rows of the
matrix; the i<span class="math">\(^{th}\)</span> process in that communicator contains the
i<span class="math">\(^{th}\)</span> block of the matrix.</p>
</section>
<section id="preallocation-of-memory-for-parallel-aij-sparse-matrices">
<h4>Preallocation of Memory for Parallel AIJ Sparse Matrices<a class="headerlink" href="#preallocation-of-memory-for-parallel-aij-sparse-matrices" title="Link to this heading">#</a></h4>
<p>As discussed above, preallocation of memory is critical for achieving
good performance during matrix assembly, as this reduces the number of
allocations and copies required. We present an example for three
processes to indicate how this may be done for the <code class="docutils notranslate"><span class="pre"><a href="../manualpages/Mat/MATMPIAIJ.html">MATMPIAIJ</a></span></code> matrix
format. Consider the 8 by 8 matrix, which is partitioned by default with
three rows on the first process, three on the second and two on the
third.</p>
<div class="math">
\[
\left( \begin{array}{cccccccccc}
1  & 2  & 0  & | & 0  & 3  & 0  & |  & 0  & 4  \\
0  & 5  & 6  & | & 7  & 0  & 0  & |  & 8  & 0 \\
9  & 0  & 10 & | & 11 & 0  & 0  & |  & 12 & 0  \\
\hline \\
13 & 0  & 14 & | & 15 & 16 & 17 & |  & 0  & 0  \\
0  & 18 & 0  & | & 19 & 20 & 21 & |  & 0  & 0 \\
0  & 0  & 0  & | & 22 & 23 & 0  & |  & 24 & 0 \\
\hline \\
25 & 26 & 27 & | & 0  & 0  & 28 & |  & 29 & 0 \\
30 & 0  & 0  & | & 31 & 32 & 33 & |  & 0  &34
\end{array} \right)
\]</div>
<p>The “diagonal” submatrix, <code class="docutils notranslate"><span class="pre">d</span></code>, on the first process is given by</p>
<div class="math">
\[
\left( \begin{array}{ccc}
1  & 2  & 0  \\
0  & 5  & 6  \\
9  & 0  & 10
\end{array} \right),
\]</div>
<p>while the “off-diagonal” submatrix, <code class="docutils notranslate"><span class="pre">o</span></code>, matrix is given by</p>
<div class="math">
\[
\left( \begin{array}{ccccc}
 0  & 3  & 0   & 0  & 4  \\
 7  & 0  & 0   & 8  & 0  \\
 11 & 0  & 0   & 12 & 0  \\
\end{array} \right).
\]</div>
<p>For the first process one could set <code class="docutils notranslate"><span class="pre">d_nz</span></code> to 2 (since each row has 2
nonzeros) or, alternatively, set <code class="docutils notranslate"><span class="pre">d_nnz</span></code> to <span class="math">\(\{2,2,2\}.\)</span> The
<code class="docutils notranslate"><span class="pre">o_nz</span></code> could be set to 2 since each row of the <code class="docutils notranslate"><span class="pre">o</span></code> matrix has 2
nonzeros, or <code class="docutils notranslate"><span class="pre">o_nnz</span></code> could be set to <span class="math">\(\{2,2,2\}\)</span>.</p>
<p>For the second process the <code class="docutils notranslate"><span class="pre">d</span></code> submatrix is given by</p>
<div class="math">
\[
\left( \begin{array}{cccccccccc}
 15 & 16 & 17 \\
 19 & 20 & 21 \\
 22 & 23 & 0
\end{array} \right) .
\]</div>
<p>Thus, one could set <code class="docutils notranslate"><span class="pre">d_nz</span></code> to 3, since the maximum number of nonzeros
in each row is 3, or alternatively one could set <code class="docutils notranslate"><span class="pre">d_nnz</span></code> to
<span class="math">\(\{3,3,2\}\)</span>, thereby indicating that the first two rows will have
3 nonzeros while the third has 2. The corresponding <code class="docutils notranslate"><span class="pre">o</span></code> submatrix for
the second process is</p>
<div class="math">
\[
\left( \begin{array}{cccccccccc}
13 & 0  & 14 &  0  & 0  \\
0  & 18 & 0  &  0  & 0 \\
0  & 0  & 0  &  24 & 0 \\
\end{array} \right)
\]</div>
<p>so that one could set <code class="docutils notranslate"><span class="pre">o_nz</span></code> to 2 or <code class="docutils notranslate"><span class="pre">o_nnz</span></code> to {2,1,1}.</p>
<p>Note that the user never directly works with the <code class="docutils notranslate"><span class="pre">d</span></code> and <code class="docutils notranslate"><span class="pre">o</span></code>
submatrices, except when preallocating storage space as indicated above.
Also, the user need not preallocate exactly the correct amount of space;
as long as a sufficiently close estimate is given, the high efficiency
for matrix assembly will remain.</p>
<p>As described above, the option <code class="docutils notranslate"><span class="pre">-info</span></code> will print information about
the success of preallocation during matrix assembly. For the
<code class="docutils notranslate"><span class="pre"><a href="../manualpages/Mat/MATMPIAIJ.html">MATMPIAIJ</a></span></code> and <code class="docutils notranslate"><span class="pre"><a href="../manualpages/Mat/MATMPIBAIJ.html">MATMPIBAIJ</a></span></code> formats, PETSc will also list the
number of elements owned by on each process that were generated on a
different process. For example, the statements</p>
<div class="highlight-c notranslate"><div class="highlight"><pre><span></span><span class="nl">MatAssemblyBegin_MPIAIJ</span><span class="p">:</span><span class="n">Stash</span><span class="w"> </span><span class="n">has</span><span class="w"> </span><span class="mi">10</span><span class="w"> </span><span class="n">entries</span><span class="p">,</span><span class="w"> </span><span class="n">uses</span><span class="w"> </span><span class="mi">0</span><span class="w"> </span><span class="n">mallocs</span>
<span class="nl">MatAssemblyBegin_MPIAIJ</span><span class="p">:</span><span class="n">Stash</span><span class="w"> </span><span class="n">has</span><span class="w"> </span><span class="mi">3</span><span class="w"> </span><span class="n">entries</span><span class="p">,</span><span class="w"> </span><span class="n">uses</span><span class="w"> </span><span class="mi">0</span><span class="w"> </span><span class="n">mallocs</span>
<span class="nl">MatAssemblyBegin_MPIAIJ</span><span class="p">:</span><span class="n">Stash</span><span class="w"> </span><span class="n">has</span><span class="w"> </span><span class="mi">5</span><span class="w"> </span><span class="n">entries</span><span class="p">,</span><span class="w"> </span><span class="n">uses</span><span class="w"> </span><span class="mi">0</span><span class="w"> </span><span class="n">mallocs</span>
</pre></div>
</div>
<p>indicate that very few values have been generated on different
processes. On the other hand, the statements</p>
<div class="highlight-c notranslate"><div class="highlight"><pre><span></span><span class="nl">MatAssemblyBegin_MPIAIJ</span><span class="p">:</span><span class="n">Stash</span><span class="w"> </span><span class="n">has</span><span class="w"> </span><span class="mi">100000</span><span class="w"> </span><span class="n">entries</span><span class="p">,</span><span class="w"> </span><span class="n">uses</span><span class="w"> </span><span class="mi">100</span><span class="w"> </span><span class="n">mallocs</span>
<span class="nl">MatAssemblyBegin_MPIAIJ</span><span class="p">:</span><span class="n">Stash</span><span class="w"> </span><span class="n">has</span><span class="w"> </span><span class="mi">77777</span><span class="w"> </span><span class="n">entries</span><span class="p">,</span><span class="w"> </span><span class="n">uses</span><span class="w"> </span><span class="mi">70</span><span class="w"> </span><span class="n">mallocs</span>
</pre></div>
</div>
<p>indicate that many values have been generated on the “wrong” processes.
This situation can be very inefficient, since the transfer of values to
the “correct” process is generally expensive. By using the command
<code class="docutils notranslate"><span class="pre"><a href="../manualpages/Mat/MatGetOwnershipRange.html">MatGetOwnershipRange</a>()</span></code> in application codes, the user should be able
to generate most entries on the owning process.</p>
<p><em>Note</em>: It is fine to generate some entries on the “wrong” process.
Often this can lead to cleaner, simpler, less buggy codes. One should
never make code overly complicated in order to generate all values
locally. Rather, one should organize the code in such a way that <em>most</em>
values are generated locally.</p>
<p>The routine <code class="docutils notranslate"><span class="pre"><a href="../manualpages/Mat/MatCreateAIJCUSPARSE.html">MatCreateAIJCUSPARSE</a>()</span></code> allows one to create GPU based matrices for NVIDIA systems.
<code class="docutils notranslate"><span class="pre"><a href="../manualpages/Mat/MatCreateAIJKokkos.html">MatCreateAIJKokkos</a>()</span></code> can create matrices for use with CPU, OpenMP, NVIDIA, AMD, or Intel based GPU systems.</p>
<p>It is sometimes difficult to compute the required preallocation information efficiently, hence PETSc provides a
special <code class="docutils notranslate"><span class="pre"><a href="../manualpages/Mat/MatType.html">MatType</a></span></code>, <code class="docutils notranslate"><span class="pre"><a href="../manualpages/Mat/MATPREALLOCATOR.html">MATPREALLOCATOR</a></span></code> that helps make computing this information more straightforward. One first creates a matrix of this type and then, using the same
code that one would use to actually compute the matrices numerical values, calls <code class="docutils notranslate"><span class="pre"><a href="../manualpages/Mat/MatSetValues.html">MatSetValues</a>()</span></code> for this matrix, without needing to provide any
preallocation information (one need not provide the matrix numerical values). Once this is complete one uses <code class="docutils notranslate"><span class="pre"><a href="../manualpages/Mat/MatPreallocatorPreallocate.html">MatPreallocatorPreallocate</a>()</span></code> to
provide the accumulated preallocation information to
the actual matrix one will use for the computations. We hope to simplify this process in the future, allowing the removal of <code class="docutils notranslate"><span class="pre"><a href="../manualpages/Mat/MATPREALLOCATOR.html">MATPREALLOCATOR</a></span></code>,
instead simply allowing the use of its efficient insertion process automatically during the first assembly of any matrix type directly without
requiring the detailed preallocation information.</p>
<p>See <a class="reference internal" href="../overview/matrix_table.html#doc-matrix"><span class="std std-ref">Summary of Matrix Types Available In PETSc</span></a> for a table of the matrix types available in PETSc.</p>
</section>
<section id="limited-memory-variable-metric-lmvm-matrices">
<span id="sec-matlmvm"></span><h4>Limited-Memory Variable Metric (LMVM) Matrices<a class="headerlink" href="#limited-memory-variable-metric-lmvm-matrices" title="Link to this heading">#</a></h4>
<p>Variable metric methods, also known as quasi-Newton methods, are
frequently used for root finding problems and approximate Jacobian
matrices or their inverses via sequential nonlinear updates based on the
secant condition. The limited-memory variants do not store the full
explicit Jacobian, and instead compute forward products and inverse
applications based on a fixed number of stored update vectors.</p>
<table class="table" id="tab-matlmvmimpl">
<caption><span class="caption-number">Table 3 </span><span class="caption-text">PETSc LMVM matrix implementations.</span><a class="headerlink" href="#tab-matlmvmimpl" title="Link to this table">#</a></caption>
<thead>
<tr class="row-odd"><th class="head"><p>Method</p></th>
<th class="head"><p>PETSc Type</p></th>
<th class="head"><p>Name</p></th>
<th class="head"><p>Property</p></th>
</tr>
</thead>
<tbody>
<tr class="row-even"><td><p>“Good” Broyden   <span id="id1">[<a class="reference internal" href="#id2239" title="Andreas Griewank. Broyden updating, the good and the bad! Optimization Stories, Documenta Mathematica. Extra Volume: Optimization Stories, pages 301–315, 2012.">ref-Gri12</a>]</span></p></td>
<td><p><code class="docutils notranslate"><span class="pre">MATLMVMBrdn</span></code></p></td>
<td><p><code class="docutils notranslate"><span class="pre">lmvmbrdn</span></code></p></td>
<td><p>Square</p></td>
</tr>
<tr class="row-odd"><td><p>“Bad” Broyden <span id="id2">[<a class="reference internal" href="#id2239" title="Andreas Griewank. Broyden updating, the good and the bad! Optimization Stories, Documenta Mathematica. Extra Volume: Optimization Stories, pages 301–315, 2012.">ref-Gri12</a>]</span></p></td>
<td><p><code class="docutils notranslate"><span class="pre">MATLMVMBadBrdn</span></code></p></td>
<td><p><code class="docutils notranslate"><span class="pre">lmvmbadbrdn</span></code></p></td>
<td><p>Square</p></td>
</tr>
<tr class="row-even"><td><p>Symmetric Rank-1 <span id="id3">[<a class="reference internal" href="#id3887" title="Jorge Nocedal and Stephen Wright. Numerical optimization. Springer Science &amp; Business Media, 2006.">ref-NW06</a>]</span></p></td>
<td><p><code class="docutils notranslate"><span class="pre">MATLMVMSR1</span></code></p></td>
<td><p><code class="docutils notranslate"><span class="pre">lmvmsr1</span></code></p></td>
<td><p>Symmetric</p></td>
</tr>
<tr class="row-odd"><td><p>Davidon-Fletcher-Powell (DFP) <span id="id4">[<a class="reference internal" href="#id3887" title="Jorge Nocedal and Stephen Wright. Numerical optimization. Springer Science &amp; Business Media, 2006.">ref-NW06</a>]</span></p></td>
<td><p><code class="docutils notranslate"><span class="pre">MATLMVMDFP</span></code></p></td>
<td><p><code class="docutils notranslate"><span class="pre">lmvmdfp</span></code></p></td>
<td><p>SPD</p></td>
</tr>
<tr class="row-even"><td><p>Dense Davidon-Fletcher-Powell (DFP) <span id="id5">[<a class="reference internal" href="#id3945" title="Jennifer B Erway and Roummel F Marcia. On solving large-scale limited-memory quasi-Newton equations. Linear Algebra and its Applications, 515:196–225, 2017.">ref-EM17b</a>]</span></p></td>
<td><p><code class="docutils notranslate"><span class="pre">MATLMVMDDFP</span></code></p></td>
<td><p><code class="docutils notranslate"><span class="pre">lmvmddfp</span></code></p></td>
<td><p>SPD</p></td>
</tr>
<tr class="row-odd"><td><p>Broyden-Fletcher-Goldfarb-Shanno (BFGS) <span id="id6">[<a class="reference internal" href="#id3887" title="Jorge Nocedal and Stephen Wright. Numerical optimization. Springer Science &amp; Business Media, 2006.">ref-NW06</a>]</span></p></td>
<td><p><code class="docutils notranslate"><span class="pre">MATLMVMBFGS</span></code></p></td>
<td><p><code class="docutils notranslate"><span class="pre">lmvmbfgs</span></code></p></td>
<td><p>SPD</p></td>
</tr>
<tr class="row-even"><td><p>Dense Broyden-Fletcher-Goldfarb-Shanno (BFGS) <span id="id7">[<a class="reference internal" href="#id3945" title="Jennifer B Erway and Roummel F Marcia. On solving large-scale limited-memory quasi-Newton equations. Linear Algebra and its Applications, 515:196–225, 2017.">ref-EM17b</a>]</span></p></td>
<td><p><code class="docutils notranslate"><span class="pre">MATLMVMDBFGS</span></code></p></td>
<td><p><code class="docutils notranslate"><span class="pre">lmvmdbfgs</span></code></p></td>
<td><p>SPD</p></td>
</tr>
<tr class="row-odd"><td><p>Dense Quasi-Newton</p></td>
<td><p><code class="docutils notranslate"><span class="pre">MATLMVMDQN</span></code></p></td>
<td><p><code class="docutils notranslate"><span class="pre">lmvmdqn</span></code></p></td>
<td><p>SPD</p></td>
</tr>
<tr class="row-even"><td><p>Restricted Broyden Family <span id="id8">[<a class="reference internal" href="#id2238" title="Jennifer B Erway and Roummel F Marcia. On solving large-scale limited-memory quasi-newton equations. Linear Algebra and its Applications, 515:196–225, 2017.">ref-EM17a</a>]</span></p></td>
<td><p><code class="docutils notranslate"><span class="pre">MATLMVMSymBrdn</span></code></p></td>
<td><p><code class="docutils notranslate"><span class="pre">lmvmsymbrdn</span></code></p></td>
<td><p>SPD</p></td>
</tr>
<tr class="row-odd"><td><p>Restricted Broyden Family (full-memory diagonal)</p></td>
<td><p><code class="docutils notranslate"><span class="pre">MATLMVMDiagBrdn</span></code></p></td>
<td><p><code class="docutils notranslate"><span class="pre">lmvmdiagbrdn</span></code></p></td>
<td><p>SPD</p></td>
</tr>
</tbody>
</table>
<p>PETSc implements seven different LMVM matrices listed in the
table above. They can be created using the
<code class="docutils notranslate"><span class="pre"><a href="../manualpages/Mat/MatCreate.html">MatCreate</a>()</span></code> and <code class="docutils notranslate"><span class="pre"><a href="../manualpages/Mat/MatSetType.html">MatSetType</a>()</span></code> workflow, and share a number of
common interface functions. We will review the most important ones
below:</p>
<ul class="simple">
<li><p><code class="docutils notranslate"><span class="pre"><a href="../manualpages/KSP/MatLMVMAllocate.html">MatLMVMAllocate</a>(<a href="../manualpages/Mat/Mat.html">Mat</a></span> <span class="pre">B,</span> <span class="pre"><a href="../manualpages/Vec/Vec.html">Vec</a></span> <span class="pre">X,</span> <span class="pre"><a href="../manualpages/Vec/Vec.html">Vec</a></span> <span class="pre">F)</span></code> – Creates the internal data
structures necessary to store nonlinear updates and compute
forward/inverse applications. The <code class="docutils notranslate"><span class="pre">X</span></code> vector defines the solution
space while the <code class="docutils notranslate"><span class="pre">F</span></code> defines the function space for the history of
updates.</p></li>
<li><p><code class="docutils notranslate"><span class="pre"><a href="../manualpages/KSP/MatLMVMUpdate.html">MatLMVMUpdate</a>(<a href="../manualpages/Mat/Mat.html">Mat</a></span> <span class="pre">B,</span> <span class="pre"><a href="../manualpages/Vec/Vec.html">Vec</a></span> <span class="pre">X,</span> <span class="pre"><a href="../manualpages/Vec/Vec.html">Vec</a></span> <span class="pre">F)</span></code> – Applies a nonlinear update to
the approximate Jacobian such that <span class="math">\(s_k = x_k - x_{k-1}\)</span> and
<span class="math">\(y_k = f(x_k) - f(x_{k-1})\)</span>, where <span class="math">\(k\)</span> is the index for
the update.</p></li>
<li><p><code class="docutils notranslate"><span class="pre"><a href="../manualpages/KSP/MatLMVMReset.html">MatLMVMReset</a>(<a href="../manualpages/Mat/Mat.html">Mat</a></span> <span class="pre">B,</span> <span class="pre"><a href="../manualpages/Sys/PetscBool.html">PetscBool</a></span> <span class="pre">destructive)</span></code> – Flushes the
accumulated nonlinear updates and resets the matrix to the initial
state. If <code class="docutils notranslate"><span class="pre">destructive</span> <span class="pre">=</span> <span class="pre"><a href="../manualpages/Sys/PETSC_TRUE.html">PETSC_TRUE</a></span></code>, the reset also destroys the
internal data structures and necessitates another allocation call
before the matrix can be updated and used for products and solves.</p></li>
<li><p><code class="docutils notranslate"><span class="pre"><a href="../manualpages/KSP/MatLMVMSetJ0.html">MatLMVMSetJ0</a>(<a href="../manualpages/Mat/Mat.html">Mat</a></span> <span class="pre">B,</span> <span class="pre"><a href="../manualpages/Mat/Mat.html">Mat</a></span> <span class="pre">J0)</span></code> – Defines the initial Jacobian to
apply the updates to. If no initial Jacobian is provided, the updates
are applied to an identity matrix.</p></li>
</ul>
<p>LMVM matrices can be applied to vectors in forward mode via
<code class="docutils notranslate"><span class="pre"><a href="../manualpages/Mat/MatMult.html">MatMult</a>()</span></code> or <code class="docutils notranslate"><span class="pre"><a href="../manualpages/Mat/MatMultAdd.html">MatMultAdd</a>()</span></code>, and in inverse mode via
<code class="docutils notranslate"><span class="pre"><a href="../manualpages/Mat/MatSolve.html">MatSolve</a>()</span></code>. They also support <code class="docutils notranslate"><span class="pre"><a href="../manualpages/Mat/MatCreateVecs.html">MatCreateVecs</a>()</span></code>, <code class="docutils notranslate"><span class="pre"><a href="../manualpages/Mat/MatDuplicate.html">MatDuplicate</a>()</span></code>
and <code class="docutils notranslate"><span class="pre"><a href="../manualpages/Mat/MatCopy.html">MatCopy</a>()</span></code> operations.</p>
<p>Restricted Broyden Family, DFP and BFGS methods, including their dense
versions, additionally implement special Jacobian initialization and
scaling options available via
<code class="docutils notranslate"><span class="pre">-mat_lmvm_scale_type</span> <span class="pre">&lt;none,scalar,diagonal&gt;</span></code>. We describe these
choices below:</p>
<ul>
<li><p><code class="docutils notranslate"><span class="pre">none</span></code> – Sets the initial Jacobian to be equal to the identity
matrix. No extra computations are required when obtaining the search
direction or updating the approximation. However, the number of
function evaluations required to converge the Newton solution is
typically much larger than what is required when using other
initializations.</p></li>
<li><p><code class="docutils notranslate"><span class="pre">scalar</span></code> – Defines the initial Jacobian as a scalar multiple of the
identity matrix. The scalar value <span class="math">\(\sigma\)</span> is chosen by solving
the one dimensional optimization problem</p>
<div class="math">
\[
  \min_\sigma \|\sigma^\alpha Y - \sigma^{\alpha - 1} S\|_F^2,
  \]</div>
<p>where <span class="math">\(S\)</span> and <span class="math">\(Y\)</span> are the matrices whose columns contain
a subset of update vectors <span class="math">\(s_k\)</span> and <span class="math">\(y_k\)</span>, and
<span class="math">\(\alpha \in [0, 1]\)</span> is defined by the user via
<code class="docutils notranslate"><span class="pre">-mat_lmvm_alpha</span></code> and has a different default value for each LMVM
implementation (e.g.: default <span class="math">\(\alpha = 1\)</span> for BFGS produces
the well-known <span class="math">\(y_k^T s_k / y_k^T y_k\)</span> scalar initialization).
The number of updates to be used in the <span class="math">\(S\)</span> and <span class="math">\(Y\)</span>
matrices is 1 by default (i.e.: the latest update only) and can be
changed via <code class="docutils notranslate"><span class="pre">-mat_lmvm_sigma_hist</span></code>. This technique is inspired by
Gilbert and Lemarechal <span id="id1">[<a class="reference internal" href="#id2905" title="J. C. Gilbert and C. Lemarechal. Some numerical experiments with variable-storage quasi-newton algorithms. Mathematical Programming, 45:407–434, 1989.">ref-GL89</a>]</span>.</p>
</li>
<li><p><code class="docutils notranslate"><span class="pre">diagonal</span></code> – Uses a full-memory restricted Broyden update formula
to construct a diagonal matrix for the Jacobian initialization.
Although the full-memory formula is utilized, the actual memory
footprint is restricted to only the vector representing the diagonal
and some additional work vectors used in its construction. The
diagonal terms are also re-scaled with every update as suggested in
<span id="id2">[<a class="reference internal" href="#id2905" title="J. C. Gilbert and C. Lemarechal. Some numerical experiments with variable-storage quasi-newton algorithms. Mathematical Programming, 45:407–434, 1989.">ref-GL89</a>]</span>. This initialization requires
the most computational effort of the available choices but typically
results in a significant reduction in the number of function
evaluations taken to compute a solution.</p></li>
</ul>
<p>The dense implementations are numerically equivalent to DFP and BFGS,
but they try to minimize memory transfer at the cost of storage
<span id="id3">[<a class="reference internal" href="#id3945" title="Jennifer B Erway and Roummel F Marcia. On solving large-scale limited-memory quasi-Newton equations. Linear Algebra and its Applications, 515:196–225, 2017.">ref-EM17b</a>]</span>. Generally, dense formulations of DFP
and BFGS, <code class="docutils notranslate"><span class="pre">MATLMVMDDFP</span></code> and <code class="docutils notranslate"><span class="pre">MATLMVMDBFGS</span></code>, should be faster than
classical recursive versions - on both CPU and GPU. It should be noted
that <code class="docutils notranslate"><span class="pre"><a href="../manualpages/Mat/MatMult.html">MatMult</a></span></code> of dense BFGS, and <code class="docutils notranslate"><span class="pre"><a href="../manualpages/Mat/MatSolve.html">MatSolve</a></span></code> of dense DFP requires
Cholesky factorization, which may be numerically unstable, if a Jacobian
option other than <code class="docutils notranslate"><span class="pre">none</span></code> is used. Therefore, the default
implementation is to enable classical recursive algorithms to avoid
the Cholesky factorization. This option can be toggled via
<code class="docutils notranslate"><span class="pre">-mat_lbfgs_recursive</span></code> and <code class="docutils notranslate"><span class="pre">-mat_ldfp_recursive</span></code>.</p>
<p>Dense Quasi-Newton, <code class="docutils notranslate"><span class="pre">MATLMVMDQN</span></code> is an implementation that uses
<code class="docutils notranslate"><span class="pre"><a href="../manualpages/Mat/MatSolve.html">MatSolve</a></span></code> of <code class="docutils notranslate"><span class="pre">MATLMVMDBFGS</span></code> for its <code class="docutils notranslate"><span class="pre"><a href="../manualpages/Mat/MatSolve.html">MatSolve</a></span></code>, and uses
<code class="docutils notranslate"><span class="pre"><a href="../manualpages/Mat/MatMult.html">MatMult</a></span></code> of <code class="docutils notranslate"><span class="pre">MATLMVMDDFP</span></code> for its <code class="docutils notranslate"><span class="pre"><a href="../manualpages/Mat/MatMult.html">MatMult</a></span></code>. It can be
seen as a hybrid implementation to avoid both recursive implementation
and Cholesky factorization, trading numerical accuracy for performances.</p>
<p>Note that the user-provided initial Jacobian via <code class="docutils notranslate"><span class="pre"><a href="../manualpages/KSP/MatLMVMSetJ0.html">MatLMVMSetJ0</a>()</span></code>
overrides and disables all built-in initialization methods.</p>
</section>
</section>
<section id="dense-matrices">
<span id="sec-matdense"></span><h3>Dense Matrices<a class="headerlink" href="#dense-matrices" title="Link to this heading">#</a></h3>
<p>PETSc provides both sequential and parallel dense matrix formats, where
each process stores its entries in a column-major array in the usual
Fortran style. To create a sequential, dense PETSc matrix, <code class="docutils notranslate"><span class="pre">A</span></code> of
dimensions <code class="docutils notranslate"><span class="pre">m</span></code> by <code class="docutils notranslate"><span class="pre">n</span></code>, the user should call</p>
<div class="highlight-c notranslate"><div class="highlight"><pre><span></span><span class="n"><a href="../manualpages/Mat/MatCreateSeqDense.html">MatCreateSeqDense</a></span><span class="p">(</span><span class="n"><a href="../manualpages/Sys/PETSC_COMM_SELF.html">PETSC_COMM_SELF</a></span><span class="p">,</span><span class="n"><a href="../manualpages/Sys/PetscInt.html">PetscInt</a></span><span class="w"> </span><span class="n">m</span><span class="p">,</span><span class="n"><a href="../manualpages/Sys/PetscInt.html">PetscInt</a></span><span class="w"> </span><span class="n">n</span><span class="p">,</span><span class="n"><a href="../manualpages/Sys/PetscScalar.html">PetscScalar</a></span><span class="w"> </span><span class="o">*</span><span class="n">data</span><span class="p">,</span><span class="n"><a href="../manualpages/Mat/Mat.html">Mat</a></span><span class="w"> </span><span class="o">*</span><span class="n">A</span><span class="p">);</span>
</pre></div>
</div>
<p>The variable <code class="docutils notranslate"><span class="pre">data</span></code> enables the user to optionally provide the
location of the data for matrix storage (intended for Fortran users who
wish to allocate their own storage space). Most users should merely set
<code class="docutils notranslate"><span class="pre">data</span></code> to <code class="docutils notranslate"><span class="pre">NULL</span></code> for PETSc to control matrix memory allocation. To
create a parallel, dense matrix, <code class="docutils notranslate"><span class="pre">A</span></code>, the user should call</p>
<div class="highlight-c notranslate"><div class="highlight"><pre><span></span><span class="n"><a href="../manualpages/Mat/MatCreateDense.html">MatCreateDense</a></span><span class="p">(</span><span class="n"><a href="../manualpages/Sys/MPI_Comm.html">MPI_Comm</a></span><span class="w"> </span><span class="n">comm</span><span class="p">,</span><span class="n"><a href="../manualpages/Sys/PetscInt.html">PetscInt</a></span><span class="w"> </span><span class="n">m</span><span class="p">,</span><span class="n"><a href="../manualpages/Sys/PetscInt.html">PetscInt</a></span><span class="w"> </span><span class="n">n</span><span class="p">,</span><span class="n"><a href="../manualpages/Sys/PetscInt.html">PetscInt</a></span><span class="w"> </span><span class="n">M</span><span class="p">,</span><span class="n"><a href="../manualpages/Sys/PetscInt.html">PetscInt</a></span><span class="w"> </span><span class="n">N</span><span class="p">,</span><span class="n"><a href="../manualpages/Sys/PetscScalar.html">PetscScalar</a></span><span class="w"> </span><span class="o">*</span><span class="n">data</span><span class="p">,</span><span class="n"><a href="../manualpages/Mat/Mat.html">Mat</a></span><span class="w"> </span><span class="o">*</span><span class="n">A</span><span class="p">)</span>
</pre></div>
</div>
<p>The arguments <code class="docutils notranslate"><span class="pre">m</span></code>, <code class="docutils notranslate"><span class="pre">n</span></code>, <code class="docutils notranslate"><span class="pre">M</span></code>, and <code class="docutils notranslate"><span class="pre">N</span></code>, indicate the number of
local rows and columns and the number of global rows and columns,
respectively. Either the local or global parameters can be replaced with
<code class="docutils notranslate"><span class="pre"><a href="../manualpages/Sys/PETSC_DECIDE.html">PETSC_DECIDE</a></span></code>, so that PETSc will determine them. The matrix is
stored with a fixed number of rows on each process, given by <code class="docutils notranslate"><span class="pre">m</span></code>, or
determined by PETSc if <code class="docutils notranslate"><span class="pre">m</span></code> is <code class="docutils notranslate"><span class="pre"><a href="../manualpages/Sys/PETSC_DECIDE.html">PETSC_DECIDE</a></span></code>.</p>
<p>PETSc does not provide parallel dense direct solvers, instead
interfacing to external packages that provide these solvers. Our focus
is on sparse iterative solvers.</p>
</section>
<section id="block-matrices">
<span id="sec-matnest"></span><h3>Block Matrices<a class="headerlink" href="#block-matrices" title="Link to this heading">#</a></h3>
<p>Block matrices arise when coupling variables with different meaning,
especially when solving problems with constraints (e.g. incompressible
flow) and “multi-physics” problems. Usually the number of blocks is
small and each block is partitioned in parallel. We illustrate for a
<span class="math">\(3\times 3\)</span> system with components labeled <span class="math">\(a,b,c\)</span>. With
some numbering of unknowns, the matrix could be written as</p>
<div class="math">
\[
\left( \begin{array}{ccc}
    A_{aa} & A_{ab} & A_{ac} \\
    A_{ba} & A_{bb} & A_{bc} \\
    A_{ca} & A_{cb} & A_{cc}
  \end{array} \right) .
\]</div>
<p>There are two fundamentally different ways that this matrix could be
stored, as a single assembled sparse matrix where entries from all
blocks are merged together (“monolithic”), or as separate assembled
matrices for each block (“nested”). These formats have different
performance characteristics depending on the operation being performed.
In particular, many preconditioners require a monolithic format, but
some that are very effective for solving block systems (see
<a class="reference internal" href="ksp.html#sec-block-matrices"><span class="std std-ref">Solving Block Matrices with PCFIELDSPLIT</span></a>) are more efficient when a nested
format is used. In order to stay flexible, we would like to be able to
use the same code to assemble block matrices in both monolithic and
nested formats. Additionally, for software maintainability and testing,
especially in a multi-physics context where different groups might be
responsible for assembling each of the blocks, it is desirable to be
able to use exactly the same code to assemble a single block
independently as to assemble it as part of a larger system. To do this,
we introduce the four spaces shown in <a class="reference internal" href="#fig-localspaces"><span class="std std-numref">Fig. 5</span></a>.</p>
<ul class="simple">
<li><p>The monolithic global space is the space in which the Krylov and
Newton solvers operate, with collective semantics across the entire
block system.</p></li>
<li><p>The split global space splits the blocks apart, but each split still
has collective semantics.</p></li>
<li><p>The split local space adds ghost points and separates the blocks.
Operations in this space can be performed with no parallel
communication. This is often the most natural, and certainly the most
powerful, space for matrix assembly code.</p></li>
<li><p>The monolithic local space can be thought of as adding ghost points
to the monolithic global space, but it is often more natural to use
it simply as a concatenation of split local spaces on each process.
It is not common to explicitly manipulate vectors or matrices in this
space (at least not during assembly), but it is a useful for
declaring which part of a matrix is being assembled.</p></li>
</ul>
<figure class="align-default" id="fig-localspaces">
<img alt="The relationship between spaces used for coupled assembly." src="../_images/localspaces.svg" /><figcaption>
<p><span class="caption-number">Fig. 5 </span><span class="caption-text">The relationship between spaces used for coupled assembly.</span><a class="headerlink" href="#fig-localspaces" title="Link to this image">#</a></p>
</figcaption>
</figure>
<p>The key to format-independent assembly is the function</p>
<div class="highlight-c notranslate"><div class="highlight"><pre><span></span><span class="n"><a href="../manualpages/Mat/MatGetLocalSubMatrix.html">MatGetLocalSubMatrix</a></span><span class="p">(</span><span class="n"><a href="../manualpages/Mat/Mat.html">Mat</a></span><span class="w"> </span><span class="n">A</span><span class="p">,</span><span class="n"><a href="../manualpages/IS/IS.html">IS</a></span><span class="w"> </span><span class="n">isrow</span><span class="p">,</span><span class="n"><a href="../manualpages/IS/IS.html">IS</a></span><span class="w"> </span><span class="n">iscol</span><span class="p">,</span><span class="n"><a href="../manualpages/Mat/Mat.html">Mat</a></span><span class="w"> </span><span class="o">*</span><span class="n">submat</span><span class="p">);</span>
</pre></div>
</div>
<p>which provides a “view” <code class="docutils notranslate"><span class="pre">submat</span></code> into a matrix <code class="docutils notranslate"><span class="pre">A</span></code> that operates in
the monolithic global space. The <code class="docutils notranslate"><span class="pre">submat</span></code> transforms from the split
local space defined by <code class="docutils notranslate"><span class="pre">iscol</span></code> to the split local space defined by
<code class="docutils notranslate"><span class="pre">isrow</span></code>. The index sets specify the parts of the monolithic local
space that <code class="docutils notranslate"><span class="pre">submat</span></code> should operate in. If a nested matrix format is
used, then <code class="docutils notranslate"><span class="pre"><a href="../manualpages/Mat/MatGetLocalSubMatrix.html">MatGetLocalSubMatrix</a>()</span></code> finds the nested block and returns
it without making any copies. In this case, <code class="docutils notranslate"><span class="pre">submat</span></code> is fully
functional and has a parallel communicator. If a monolithic matrix
format is used, then <code class="docutils notranslate"><span class="pre"><a href="../manualpages/Mat/MatGetLocalSubMatrix.html">MatGetLocalSubMatrix</a>()</span></code> returns a proxy matrix
on <code class="docutils notranslate"><span class="pre"><a href="../manualpages/Sys/PETSC_COMM_SELF.html">PETSC_COMM_SELF</a></span></code> that does not provide values or implement
<code class="docutils notranslate"><span class="pre"><a href="../manualpages/Mat/MatMult.html">MatMult</a>()</span></code>, but does implement <code class="docutils notranslate"><span class="pre"><a href="../manualpages/Mat/MatSetValuesLocal.html">MatSetValuesLocal</a>()</span></code> and, if
<code class="docutils notranslate"><span class="pre">isrow,iscol</span></code> have a constant block size,
<code class="docutils notranslate"><span class="pre"><a href="../manualpages/Mat/MatSetValuesBlockedLocal.html">MatSetValuesBlockedLocal</a>()</span></code>. Note that although <code class="docutils notranslate"><span class="pre">submat</span></code> may not be
a fully functional matrix and the caller does not even know a priori
which communicator it will reside on, it always implements the local
assembly functions (which are not collective). The index sets
<code class="docutils notranslate"><span class="pre">isrow,iscol</span></code> can be obtained using <code class="docutils notranslate"><span class="pre"><a href="../manualpages/DMComposite/DMCompositeGetLocalISs.html">DMCompositeGetLocalISs</a>()</span></code> if
<code class="docutils notranslate"><span class="pre"><a href="../manualpages/DMComposite/DMCOMPOSITE.html">DMCOMPOSITE</a></span></code> is being used. <code class="docutils notranslate"><span class="pre"><a href="../manualpages/DMComposite/DMCOMPOSITE.html">DMCOMPOSITE</a></span></code> can also be used to create
matrices, in which case the <code class="docutils notranslate"><span class="pre"><a href="../manualpages/Mat/MATNEST.html">MATNEST</a></span></code> format can be specified using
<code class="docutils notranslate"><span class="pre">-prefix_dm_mat_type</span> <span class="pre">nest</span></code> and <code class="docutils notranslate"><span class="pre"><a href="../manualpages/Mat/MATAIJ.html">MATAIJ</a></span></code> can be specified using
<code class="docutils notranslate"><span class="pre">-prefix_dm_mat_type</span> <span class="pre">aij</span></code>. See
<a href="../src/snes/tutorials/ex28.c.html">SNES Tutorial ex28</a>
for a simple example using this interface.</p>
</section>
</section>
<section id="basic-matrix-operations">
<span id="sec-matoptions"></span><h2>Basic Matrix Operations<a class="headerlink" href="#basic-matrix-operations" title="Link to this heading">#</a></h2>
<p>Table <a class="reference internal" href="#fig-matrixops"><span class="std std-ref">2.2</span></a> summarizes basic PETSc matrix operations.
We briefly discuss a few of these routines in more detail below.</p>
<p>The parallel matrix can multiply a vector with <code class="docutils notranslate"><span class="pre">n</span></code> local entries,
returning a vector with <code class="docutils notranslate"><span class="pre">m</span></code> local entries. That is, to form the
product</p>
<div class="highlight-c notranslate"><div class="highlight"><pre><span></span><span class="n"><a href="../manualpages/Mat/MatMult.html">MatMult</a></span><span class="p">(</span><span class="n"><a href="../manualpages/Mat/Mat.html">Mat</a></span><span class="w"> </span><span class="n">A</span><span class="p">,</span><span class="n"><a href="../manualpages/Vec/Vec.html">Vec</a></span><span class="w"> </span><span class="n">x</span><span class="p">,</span><span class="n"><a href="../manualpages/Vec/Vec.html">Vec</a></span><span class="w"> </span><span class="n">y</span><span class="p">);</span>
</pre></div>
</div>
<p>the vectors <code class="docutils notranslate"><span class="pre">x</span></code> and <code class="docutils notranslate"><span class="pre">y</span></code> should be generated with</p>
<div class="highlight-c notranslate"><div class="highlight"><pre><span></span><span class="n"><a href="../manualpages/Vec/VecCreateMPI.html">VecCreateMPI</a></span><span class="p">(</span><span class="n"><a href="../manualpages/Sys/MPI_Comm.html">MPI_Comm</a></span><span class="w"> </span><span class="n">comm</span><span class="p">,</span><span class="n">n</span><span class="p">,</span><span class="n">N</span><span class="p">,</span><span class="o">&amp;</span><span class="n">x</span><span class="p">);</span>
<span class="n"><a href="../manualpages/Vec/VecCreateMPI.html">VecCreateMPI</a></span><span class="p">(</span><span class="n"><a href="../manualpages/Sys/MPI_Comm.html">MPI_Comm</a></span><span class="w"> </span><span class="n">comm</span><span class="p">,</span><span class="n">m</span><span class="p">,</span><span class="n">M</span><span class="p">,</span><span class="o">&amp;</span><span class="n">y</span><span class="p">);</span>
</pre></div>
</div>
<p>By default, if the user lets PETSc decide the number of components to be
stored locally (by passing in <code class="docutils notranslate"><span class="pre"><a href="../manualpages/Sys/PETSC_DECIDE.html">PETSC_DECIDE</a></span></code> as the second argument to
<code class="docutils notranslate"><span class="pre"><a href="../manualpages/Vec/VecCreateMPI.html">VecCreateMPI</a>()</span></code> or using <code class="docutils notranslate"><span class="pre"><a href="../manualpages/Vec/VecCreate.html">VecCreate</a>()</span></code>), vectors and matrices of
the same dimension are automatically compatible for parallel
matrix-vector operations.</p>
<p>Along with the matrix-vector multiplication routine, there is a version
for the transpose of the matrix,</p>
<div class="highlight-c notranslate"><div class="highlight"><pre><span></span><span class="n"><a href="../manualpages/Mat/MatMultTranspose.html">MatMultTranspose</a></span><span class="p">(</span><span class="n"><a href="../manualpages/Mat/Mat.html">Mat</a></span><span class="w"> </span><span class="n">A</span><span class="p">,</span><span class="n"><a href="../manualpages/Vec/Vec.html">Vec</a></span><span class="w"> </span><span class="n">x</span><span class="p">,</span><span class="n"><a href="../manualpages/Vec/Vec.html">Vec</a></span><span class="w"> </span><span class="n">y</span><span class="p">);</span>
</pre></div>
</div>
<p>There are also versions that add the result to another vector:</p>
<div class="highlight-c notranslate"><div class="highlight"><pre><span></span><span class="n"><a href="../manualpages/Mat/MatMultAdd.html">MatMultAdd</a></span><span class="p">(</span><span class="n"><a href="../manualpages/Mat/Mat.html">Mat</a></span><span class="w"> </span><span class="n">A</span><span class="p">,</span><span class="n"><a href="../manualpages/Vec/Vec.html">Vec</a></span><span class="w"> </span><span class="n">x</span><span class="p">,</span><span class="n"><a href="../manualpages/Vec/Vec.html">Vec</a></span><span class="w"> </span><span class="n">y</span><span class="p">,</span><span class="n"><a href="../manualpages/Vec/Vec.html">Vec</a></span><span class="w"> </span><span class="n">w</span><span class="p">);</span>
<span class="n"><a href="../manualpages/Mat/MatMultTransposeAdd.html">MatMultTransposeAdd</a></span><span class="p">(</span><span class="n"><a href="../manualpages/Mat/Mat.html">Mat</a></span><span class="w"> </span><span class="n">A</span><span class="p">,</span><span class="n"><a href="../manualpages/Vec/Vec.html">Vec</a></span><span class="w"> </span><span class="n">x</span><span class="p">,</span><span class="n"><a href="../manualpages/Vec/Vec.html">Vec</a></span><span class="w"> </span><span class="n">y</span><span class="p">,</span><span class="n"><a href="../manualpages/Vec/Vec.html">Vec</a></span><span class="w"> </span><span class="n">w</span><span class="p">);</span>
</pre></div>
</div>
<p>These routines, respectively, produce <span class="math">\(w = A*x + y\)</span> and
<span class="math">\(w = A^{T}*x + y\)</span> . In C it is legal for the vectors <code class="docutils notranslate"><span class="pre">y</span></code> and
<code class="docutils notranslate"><span class="pre">w</span></code> to be identical. In Fortran, this situation is forbidden by the
language standard, but we allow it anyway.</p>
<p>One can print a matrix (sequential or parallel) to the screen with the
command</p>
<div class="highlight-c notranslate"><div class="highlight"><pre><span></span><span class="n"><a href="../manualpages/Mat/MatView.html">MatView</a></span><span class="p">(</span><span class="n"><a href="../manualpages/Mat/Mat.html">Mat</a></span><span class="w"> </span><span class="n">mat</span><span class="p">,</span><span class="n"><a href="../manualpages/Viewer/PETSC_VIEWER_STDOUT_WORLD.html">PETSC_VIEWER_STDOUT_WORLD</a></span><span class="p">);</span>
</pre></div>
</div>
<p>Other viewers can be used as well. For instance, one can draw the
nonzero structure of the matrix into the default X-window with the
command</p>
<div class="highlight-c notranslate"><div class="highlight"><pre><span></span><span class="n"><a href="../manualpages/Mat/MatView.html">MatView</a></span><span class="p">(</span><span class="n"><a href="../manualpages/Mat/Mat.html">Mat</a></span><span class="w"> </span><span class="n">mat</span><span class="p">,</span><span class="n"><a href="../manualpages/Viewer/PETSC_VIEWER_DRAW_WORLD.html">PETSC_VIEWER_DRAW_WORLD</a></span><span class="p">);</span>
</pre></div>
</div>
<p>Also one can use</p>
<div class="highlight-c notranslate"><div class="highlight"><pre><span></span><span class="n"><a href="../manualpages/Mat/MatView.html">MatView</a></span><span class="p">(</span><span class="n"><a href="../manualpages/Mat/Mat.html">Mat</a></span><span class="w"> </span><span class="n">mat</span><span class="p">,</span><span class="n"><a href="../manualpages/Viewer/PetscViewer.html">PetscViewer</a></span><span class="w"> </span><span class="n">viewer</span><span class="p">);</span>
</pre></div>
</div>
<p>where <code class="docutils notranslate"><span class="pre">viewer</span></code> was obtained with <code class="docutils notranslate"><span class="pre"><a href="../manualpages/Viewer/PetscViewerDrawOpen.html">PetscViewerDrawOpen</a>()</span></code>. Additional
viewers and options are given in the <code class="docutils notranslate"><span class="pre"><a href="../manualpages/Mat/MatView.html">MatView</a>()</span></code> man page and
<a class="reference internal" href="other.html#sec-viewers"><span class="std std-ref">Viewers: Looking at PETSc Objects</span></a>.</p>
<table class="table" id="fig-matrixops">
<caption><span class="caption-number">Table 4 </span><span class="caption-text">PETSc Matrix Operations</span><a class="headerlink" href="#fig-matrixops" title="Link to this table">#</a></caption>
<thead>
<tr class="row-odd"><th class="head"><p>Function Name</p></th>
<th class="head"><p>Operation</p></th>
</tr>
</thead>
<tbody>
<tr class="row-even"><td><p><code class="docutils notranslate"><span class="pre"><a href="../manualpages/Mat/MatAXPY.html">MatAXPY</a>(<a href="../manualpages/Mat/Mat.html">Mat</a></span> <span class="pre">Y,</span> <span class="pre"><a href="../manualpages/Sys/PetscScalar.html">PetscScalar</a></span> <span class="pre">a,</span> <span class="pre"><a href="../manualpages/Mat/Mat.html">Mat</a></span> <span class="pre">X,</span> <span class="pre"><a href="../manualpages/Mat/MatStructure.html">MatStructure</a></span> <span class="pre">s);</span></code></p></td>
<td><p><span class="math">\(Y = Y + a*X\)</span></p></td>
</tr>
<tr class="row-odd"><td><p><code class="docutils notranslate"><span class="pre"><a href="../manualpages/Mat/MatAYPX.html">MatAYPX</a>(<a href="../manualpages/Mat/Mat.html">Mat</a></span> <span class="pre">Y,</span> <span class="pre"><a href="../manualpages/Sys/PetscScalar.html">PetscScalar</a></span> <span class="pre">a,</span> <span class="pre"><a href="../manualpages/Mat/Mat.html">Mat</a></span> <span class="pre">X,</span> <span class="pre"><a href="../manualpages/Mat/MatStructure.html">MatStructure</a></span> <span class="pre">s);</span></code></p></td>
<td><p><span class="math">\(Y = a*Y + X\)</span></p></td>
</tr>
<tr class="row-even"><td><p><code class="docutils notranslate"><span class="pre"><a href="../manualpages/Mat/MatMult.html">MatMult</a>(<a href="../manualpages/Mat/Mat.html">Mat</a></span> <span class="pre">A,<a href="../manualpages/Vec/Vec.html">Vec</a></span> <span class="pre">x,</span> <span class="pre"><a href="../manualpages/Vec/Vec.html">Vec</a></span> <span class="pre">y);</span></code></p></td>
<td><p><span class="math">\(y = A*x\)</span></p></td>
</tr>
<tr class="row-odd"><td><p><code class="docutils notranslate"><span class="pre"><a href="../manualpages/Mat/MatMultAdd.html">MatMultAdd</a>(<a href="../manualpages/Mat/Mat.html">Mat</a></span> <span class="pre">A,<a href="../manualpages/Vec/Vec.html">Vec</a></span> <span class="pre">x,</span> <span class="pre"><a href="../manualpages/Vec/Vec.html">Vec</a></span> <span class="pre">y,<a href="../manualpages/Vec/Vec.html">Vec</a></span> <span class="pre">z);</span></code></p></td>
<td><p><span class="math">\(z = y + A*x\)</span></p></td>
</tr>
<tr class="row-even"><td><p><code class="docutils notranslate"><span class="pre"><a href="../manualpages/Mat/MatMultTranspose.html">MatMultTranspose</a>(<a href="../manualpages/Mat/Mat.html">Mat</a></span> <span class="pre">A,<a href="../manualpages/Vec/Vec.html">Vec</a></span> <span class="pre">x,</span> <span class="pre"><a href="../manualpages/Vec/Vec.html">Vec</a></span> <span class="pre">y);</span></code></p></td>
<td><p><span class="math">\(y = A^{T}*x\)</span></p></td>
</tr>
<tr class="row-odd"><td><p><code class="docutils notranslate"><span class="pre"><a href="../manualpages/Mat/MatMultTransposeAdd.html">MatMultTransposeAdd</a>(<a href="../manualpages/Mat/Mat.html">Mat</a></span> <span class="pre">A,</span> <span class="pre"><a href="../manualpages/Vec/Vec.html">Vec</a></span> <span class="pre">x,</span> <span class="pre"><a href="../manualpages/Vec/Vec.html">Vec</a></span> <span class="pre">y,</span> <span class="pre"><a href="../manualpages/Vec/Vec.html">Vec</a></span> <span class="pre">z);</span></code></p></td>
<td><p><span class="math">\(z = y + A^{T}*x\)</span></p></td>
</tr>
<tr class="row-even"><td><p><code class="docutils notranslate"><span class="pre"><a href="../manualpages/Mat/MatNorm.html">MatNorm</a>(<a href="../manualpages/Mat/Mat.html">Mat</a></span> <span class="pre">A,<a href="../manualpages/Vec/NormType.html">NormType</a></span> <span class="pre">type,</span> <span class="pre"><a href="../manualpages/Sys/PetscReal.html">PetscReal</a></span> <span class="pre">*r);</span></code></p></td>
<td><p><span class="math">\(r = A_{type}\)</span></p></td>
</tr>
<tr class="row-odd"><td><p><code class="docutils notranslate"><span class="pre"><a href="../manualpages/Mat/MatDiagonalScale.html">MatDiagonalScale</a>(<a href="../manualpages/Mat/Mat.html">Mat</a></span> <span class="pre">A,<a href="../manualpages/Vec/Vec.html">Vec</a></span> <span class="pre">l,<a href="../manualpages/Vec/Vec.html">Vec</a></span> <span class="pre">r);</span></code></p></td>
<td><p><span class="math">\(A = \text{diag}(l)*A*\text{diag}(r)\)</span></p></td>
</tr>
<tr class="row-even"><td><p><code class="docutils notranslate"><span class="pre"><a href="../manualpages/Mat/MatScale.html">MatScale</a>(<a href="../manualpages/Mat/Mat.html">Mat</a></span> <span class="pre">A,<a href="../manualpages/Sys/PetscScalar.html">PetscScalar</a></span> <span class="pre">a);</span></code></p></td>
<td><p><span class="math">\(A = a*A\)</span></p></td>
</tr>
<tr class="row-odd"><td><p><code class="docutils notranslate"><span class="pre"><a href="../manualpages/Mat/MatConvert.html">MatConvert</a>(<a href="../manualpages/Mat/Mat.html">Mat</a></span> <span class="pre">A,</span> <span class="pre"><a href="../manualpages/Mat/MatType.html">MatType</a></span> <span class="pre">type,</span> <span class="pre"><a href="../manualpages/Mat/Mat.html">Mat</a></span> <span class="pre">*B);</span></code></p></td>
<td><p><span class="math">\(B = A\)</span></p></td>
</tr>
<tr class="row-even"><td><p><code class="docutils notranslate"><span class="pre"><a href="../manualpages/Mat/MatCopy.html">MatCopy</a>(<a href="../manualpages/Mat/Mat.html">Mat</a></span> <span class="pre">A,</span> <span class="pre"><a href="../manualpages/Mat/Mat.html">Mat</a></span> <span class="pre">B,</span> <span class="pre"><a href="../manualpages/Mat/MatStructure.html">MatStructure</a></span> <span class="pre">s);</span></code></p></td>
<td><p><span class="math">\(B = A\)</span></p></td>
</tr>
<tr class="row-odd"><td><p><code class="docutils notranslate"><span class="pre"><a href="../manualpages/Mat/MatGetDiagonal.html">MatGetDiagonal</a>(<a href="../manualpages/Mat/Mat.html">Mat</a></span> <span class="pre">A,</span> <span class="pre"><a href="../manualpages/Vec/Vec.html">Vec</a></span> <span class="pre">x);</span></code></p></td>
<td><p><span class="math">\(x = \text{diag}(A)\)</span></p></td>
</tr>
<tr class="row-even"><td><p><code class="docutils notranslate"><span class="pre"><a href="../manualpages/Mat/MatTranspose.html">MatTranspose</a>(<a href="../manualpages/Mat/Mat.html">Mat</a></span> <span class="pre">A,</span> <span class="pre"><a href="../manualpages/Mat/MatReuse.html">MatReuse</a>,</span> <span class="pre"><a href="../manualpages/Mat/Mat.html">Mat</a>*</span> <span class="pre">B);</span></code></p></td>
<td><p><span class="math">\(B = A^{T}\)</span></p></td>
</tr>
<tr class="row-odd"><td><p><code class="docutils notranslate"><span class="pre"><a href="../manualpages/Mat/MatZeroEntries.html">MatZeroEntries</a>(<a href="../manualpages/Mat/Mat.html">Mat</a></span> <span class="pre">A);</span></code></p></td>
<td><p><span class="math">\(A = 0\)</span></p></td>
</tr>
<tr class="row-even"><td><p><code class="docutils notranslate"><span class="pre"><a href="../manualpages/Mat/MatShift.html">MatShift</a>(<a href="../manualpages/Mat/Mat.html">Mat</a></span> <span class="pre">Y,</span> <span class="pre"><a href="../manualpages/Sys/PetscScalar.html">PetscScalar</a></span> <span class="pre">a);</span></code></p></td>
<td><p><span class="math">\(Y =  Y + a*I\)</span></p></td>
</tr>
</tbody>
</table>
<table class="table" id="fig-matstructure">
<caption><span class="caption-number">Table 5 </span><span class="caption-text">Values of MatStructure</span><a class="headerlink" href="#fig-matstructure" title="Link to this table">#</a></caption>
<thead>
<tr class="row-odd"><th class="head"><p>Name</p></th>
<th class="head"><p>Meaning</p></th>
</tr>
</thead>
<tbody>
<tr class="row-even"><td><p><code class="docutils notranslate"><span class="pre"><a href="../manualpages/Mat/MatStructure.html">SAME_NONZERO_PATTERN</a></span></code></p></td>
<td><p>the matrices have an identical nonzero pattern</p></td>
</tr>
<tr class="row-odd"><td><p><code class="docutils notranslate"><span class="pre"><a href="../manualpages/Mat/MatStructure.html">DIFFERENT_NONZERO_PATTERN</a></span></code></p></td>
<td><p>the matrices may have a different nonzero pattern</p></td>
</tr>
<tr class="row-even"><td><p><code class="docutils notranslate"><span class="pre"><a href="../manualpages/Mat/MatStructure.html">SUBSET_NONZERO_PATTERN</a></span></code></p></td>
<td><p>the second matrix has a subset of the nonzeros in the first matrix</p></td>
</tr>
<tr class="row-odd"><td><p><code class="docutils notranslate"><span class="pre"><a href="../manualpages/Mat/MatStructure.html">UNKNOWN_NONZERO_PATTERN</a></span></code></p></td>
<td><p>there is nothing known about the relation between the nonzero patterns of the two matrices</p></td>
</tr>
</tbody>
</table>
<p>The <code class="docutils notranslate"><span class="pre"><a href="../manualpages/Vec/NormType.html">NormType</a></span></code> argument to <code class="docutils notranslate"><span class="pre"><a href="../manualpages/Mat/MatNorm.html">MatNorm</a>()</span></code> is one of <code class="docutils notranslate"><span class="pre"><a href="../manualpages/Vec/NORM_1.html">NORM_1</a></span></code>,
<code class="docutils notranslate"><span class="pre"><a href="../manualpages/Vec/NORM_INFINITY.html">NORM_INFINITY</a></span></code>, and <code class="docutils notranslate"><span class="pre"><a href="../manualpages/Vec/NORM_FROBENIUS.html">NORM_FROBENIUS</a></span></code>.</p>
</section>
<section id="application-specific-custom-matrices">
<span id="sec-matrixfree"></span><h2>Application Specific Custom Matrices<a class="headerlink" href="#application-specific-custom-matrices" title="Link to this heading">#</a></h2>
<p>Some people like to use matrix-free methods, which do
not require explicit storage of the matrix, for the numerical solution
of partial differential equations.
Similarly, users may already have a custom matrix data structure and routines
for that data structure and would like to wrap their code up into a <code class="docutils notranslate"><span class="pre"><a href="../manualpages/Mat/Mat.html">Mat</a></span></code>;
that is, provide their own custom matrix type.</p>
<p>To use the PETSc provided matrix-free matrix that uses finite differencing to approximate the matrix-vector product
use <code class="docutils notranslate"><span class="pre"><a href="../manualpages/Mat/MatCreateMFFD.html">MatCreateMFFD</a>()</span></code>, see <a class="reference internal" href="snes.html#sec-nlmatrixfree"><span class="std std-ref">Matrix-Free Methods</span></a>.
To provide your own matrix operations (such as <code class="docutils notranslate"><span class="pre"><a href="../manualpages/Mat/MatMult.html">MatMult</a>()</span></code>)
use the following command to create a <code class="docutils notranslate"><span class="pre"><a href="../manualpages/Mat/Mat.html">Mat</a></span></code> structure
without ever actually generating the matrix:</p>
<div class="highlight-c notranslate"><div class="highlight"><pre><span></span><span class="n"><a href="../manualpages/Mat/MatCreateShell.html">MatCreateShell</a></span><span class="p">(</span><span class="n"><a href="../manualpages/Sys/MPI_Comm.html">MPI_Comm</a></span><span class="w"> </span><span class="n">comm</span><span class="p">,</span><span class="n"><a href="../manualpages/Sys/PetscInt.html">PetscInt</a></span><span class="w"> </span><span class="n">m</span><span class="p">,</span><span class="n"><a href="../manualpages/Sys/PetscInt.html">PetscInt</a></span><span class="w"> </span><span class="n">n</span><span class="p">,</span><span class="n"><a href="../manualpages/Sys/PetscInt.html">PetscInt</a></span><span class="w"> </span><span class="n">M</span><span class="p">,</span><span class="n"><a href="../manualpages/Sys/PetscInt.html">PetscInt</a></span><span class="w"> </span><span class="n">N</span><span class="p">,</span><span class="kt">void</span><span class="w"> </span><span class="o">*</span><span class="n">ctx</span><span class="p">,</span><span class="n"><a href="../manualpages/Mat/Mat.html">Mat</a></span><span class="w"> </span><span class="o">*</span><span class="n">mat</span><span class="p">);</span>
</pre></div>
</div>
<p>Here <code class="docutils notranslate"><span class="pre">M</span></code> and <code class="docutils notranslate"><span class="pre">N</span></code> are the global matrix dimensions (rows and
columns), <code class="docutils notranslate"><span class="pre">m</span></code> and <code class="docutils notranslate"><span class="pre">n</span></code> are the local matrix dimensions, and <code class="docutils notranslate"><span class="pre">ctx</span></code>
is a pointer to data needed by any user-defined shell matrix operations;
the manual page has additional details about these parameters. Most
matrix-free algorithms require only the application of the linear
operator to a vector. To provide this action, the user must write a
routine with the calling sequence</p>
<div class="highlight-c notranslate"><div class="highlight"><pre><span></span><span class="n">UserMult</span><span class="p">(</span><span class="n"><a href="../manualpages/Mat/Mat.html">Mat</a></span><span class="w"> </span><span class="n">mat</span><span class="p">,</span><span class="n"><a href="../manualpages/Vec/Vec.html">Vec</a></span><span class="w"> </span><span class="n">x</span><span class="p">,</span><span class="n"><a href="../manualpages/Vec/Vec.html">Vec</a></span><span class="w"> </span><span class="n">y</span><span class="p">);</span>
</pre></div>
</div>
<p>and then associate it with the matrix, <code class="docutils notranslate"><span class="pre">mat</span></code>, by using the command</p>
<div class="highlight-c notranslate"><div class="highlight"><pre><span></span><span class="n"><a href="../manualpages/Mat/MatShellSetOperation.html">MatShellSetOperation</a></span><span class="p">(</span><span class="n"><a href="../manualpages/Mat/Mat.html">Mat</a></span><span class="w"> </span><span class="n">mat</span><span class="p">,</span><span class="n">MatOperation</span><span class="w"> </span><span class="n">MATOP_MULT</span><span class="p">,</span><span class="w"> </span><span class="p">(</span><span class="kt">void</span><span class="p">(</span><span class="o">*</span><span class="p">)(</span><span class="kt">void</span><span class="p">))</span><span class="w"> </span><span class="n"><a href="../manualpages/Sys/PetscErrorCode.html">PetscErrorCode</a></span><span class="w"> </span><span class="p">(</span><span class="o">*</span><span class="n">UserMult</span><span class="p">)(</span><span class="n"><a href="../manualpages/Mat/Mat.html">Mat</a></span><span class="p">,</span><span class="n"><a href="../manualpages/Vec/Vec.html">Vec</a></span><span class="p">,</span><span class="n"><a href="../manualpages/Vec/Vec.html">Vec</a></span><span class="p">));</span>
</pre></div>
</div>
<p>Here <code class="docutils notranslate"><span class="pre">MATOP_MULT</span></code> is the name of the operation for matrix-vector
multiplication. Within each user-defined routine (such as
<code class="docutils notranslate"><span class="pre">UserMult()</span></code>), the user should call <code class="docutils notranslate"><span class="pre"><a href="../manualpages/Mat/MatShellGetContext.html">MatShellGetContext</a>()</span></code> to obtain
the user-defined context, <code class="docutils notranslate"><span class="pre">ctx</span></code>, that was set by <code class="docutils notranslate"><span class="pre"><a href="../manualpages/Mat/MatCreateShell.html">MatCreateShell</a>()</span></code>.
This shell matrix can be used with the iterative linear equation solvers
discussed in the following chapters.</p>
<p>The routine <code class="docutils notranslate"><span class="pre"><a href="../manualpages/Mat/MatShellSetOperation.html">MatShellSetOperation</a>()</span></code> can be used to set any other
matrix operations as well. The file
<code class="docutils notranslate"><span class="pre">$PETSC_DIR/include/petscmat.h</span></code> (<a href="../include/petscmat.h.html">source</a>)
provides a complete list of matrix operations, which have the form
<code class="docutils notranslate"><span class="pre">MATOP_&lt;OPERATION&gt;</span></code>, where <code class="docutils notranslate"><span class="pre">&lt;OPERATION&gt;</span></code> is the name (in all capital
letters) of the user interface routine (for example, <code class="docutils notranslate"><span class="pre"><a href="../manualpages/Mat/MatMult.html">MatMult</a>()</span></code>
<span class="math">\(\to\)</span> <code class="docutils notranslate"><span class="pre">MATOP_MULT</span></code>). All user-provided functions have the same
calling sequence as the usual matrix interface routines, since the
user-defined functions are intended to be accessed through the same
interface, e.g., <code class="docutils notranslate"><span class="pre"><a href="../manualpages/Mat/MatMult.html">MatMult</a>(<a href="../manualpages/Mat/Mat.html">Mat</a>,<a href="../manualpages/Vec/Vec.html">Vec</a>,<a href="../manualpages/Vec/Vec.html">Vec</a>)</span></code> <span class="math">\(\to\)</span>
<code class="docutils notranslate"><span class="pre">UserMult(<a href="../manualpages/Mat/Mat.html">Mat</a>,<a href="../manualpages/Vec/Vec.html">Vec</a>,<a href="../manualpages/Vec/Vec.html">Vec</a>)</span></code>. The final argument for
<code class="docutils notranslate"><span class="pre"><a href="../manualpages/Mat/MatShellSetOperation.html">MatShellSetOperation</a>()</span></code> needs to be cast to a <code class="docutils notranslate"><span class="pre">void</span> <span class="pre">*</span></code>, since the
final argument could (depending on the <code class="docutils notranslate"><span class="pre">MatOperation</span></code>) be a variety of
different functions.</p>
<p>Note that <code class="docutils notranslate"><span class="pre"><a href="../manualpages/Mat/MatShellSetOperation.html">MatShellSetOperation</a>()</span></code> can also be used as a “backdoor”
means of introducing user-defined changes in matrix operations for other
storage formats (for example, to override the default LU factorization
routine supplied within PETSc for the <code class="docutils notranslate"><span class="pre"><a href="../manualpages/Mat/MATSEQAIJ.html">MATSEQAIJ</a></span></code> format). However, we
urge anyone who introduces such changes to use caution, since it would
be very easy to accidentally create a bug in the new routine that could
affect other routines as well.</p>
<p>See also <a class="reference internal" href="snes.html#sec-nlmatrixfree"><span class="std std-ref">Matrix-Free Methods</span></a> for details on one set of
helpful utilities for using the matrix-free approach for nonlinear
solvers.</p>
</section>
<section id="transposes-of-matrices">
<span id="sec-mattranspose"></span><h2>Transposes of Matrices<a class="headerlink" href="#transposes-of-matrices" title="Link to this heading">#</a></h2>
<p>PETSc provides several ways to work with transposes of matrix.</p>
<div class="highlight-c notranslate"><div class="highlight"><pre><span></span><span class="n"><a href="../manualpages/Mat/MatTranspose.html">MatTranspose</a></span><span class="p">(</span><span class="n"><a href="../manualpages/Mat/Mat.html">Mat</a></span><span class="w"> </span><span class="n">A</span><span class="p">,</span><span class="n"><a href="../manualpages/Mat/MatReuse.html">MatReuse</a></span><span class="w"> </span><span class="n"><a href="../manualpages/Mat/MatReuse.html">MAT_INITIAL_MATRIX</a></span><span class="w"> </span><span class="n">or</span><span class="w"> </span><span class="n"><a href="../manualpages/Mat/MatReuse.html">MAT_INPLACE_MATRIX</a></span><span class="w"> </span><span class="n">or</span><span class="w"> </span><span class="n"><a href="../manualpages/Mat/MatReuse.html">MAT_REUSE_MATRIX</a></span><span class="p">,</span><span class="n"><a href="../manualpages/Mat/Mat.html">Mat</a></span><span class="w"> </span><span class="o">*</span><span class="n">B</span><span class="p">)</span>
</pre></div>
</div>
<p>will either do an in-place or out-of-place matrix explicit formation of the matrix transpose. After it has been called
with <code class="docutils notranslate"><span class="pre"><a href="../manualpages/Mat/MatReuse.html">MAT_INPLACE_MATRIX</a></span></code> it may be called again with <code class="docutils notranslate"><span class="pre"><a href="../manualpages/Mat/MatReuse.html">MAT_REUSE_MATRIX</a></span></code> and it will recompute the transpose if the A
matrix has changed. Internally it keeps track of whether the nonzero pattern of A has not changed so
will reuse the symbolic transpose when possible for efficiency.</p>
<div class="highlight-c notranslate"><div class="highlight"><pre><span></span><span class="n"><a href="../manualpages/Mat/MatTransposeSymbolic.html">MatTransposeSymbolic</a></span><span class="p">(</span><span class="n"><a href="../manualpages/Mat/Mat.html">Mat</a></span><span class="w"> </span><span class="n">A</span><span class="p">,</span><span class="n"><a href="../manualpages/Mat/Mat.html">Mat</a></span><span class="w"> </span><span class="o">*</span><span class="n">B</span><span class="p">)</span>
</pre></div>
</div>
<p>only does the symbolic transpose on the matrix. After it is called <code class="docutils notranslate"><span class="pre"><a href="../manualpages/Mat/MatTranspose.html">MatTranspose</a>()</span></code> may be called with
<code class="docutils notranslate"><span class="pre"><a href="../manualpages/Mat/MatReuse.html">MAT_REUSE_MATRIX</a></span></code> to compute the numerical transpose.</p>
<p>Occasionally one may already have a B matrix with the needed sparsity pattern to store the transpose and wants to reuse that
space instead of creating a new matrix by calling <code class="docutils notranslate"><span class="pre"><a href="../manualpages/Mat/MatTranspose.html">MatTranspose</a></span></code>(A,``MAT_INITIAL_MATRIX``,&amp;B) but they cannot just call
<code class="docutils notranslate"><span class="pre"><a href="../manualpages/Mat/MatTranspose.html">MatTranspose</a></span></code>(A,``MAT_REUSE_MATRIX``,&amp;B) so instead they can call <code class="docutils notranslate"><span class="pre">MatTransposeSetPrecusor</span></code>(A,B) and then call
<code class="docutils notranslate"><span class="pre"><a href="../manualpages/Mat/MatTranspose.html">MatTranspose</a></span></code>(A,``MAT_REUSE_MATRIX``,&amp;B). This routine just provides to B the meta-data it needs to compute the numerical
factorization efficiently.</p>
<p>The routine <code class="docutils notranslate"><span class="pre"><a href="../manualpages/Mat/MatCreateTranspose.html">MatCreateTranspose</a></span></code>(A,&amp;B) provides a surrogate matrix B that behaviors like the transpose of A without forming
the transpose explicitly. For example, <code class="docutils notranslate"><span class="pre"><a href="../manualpages/Mat/MatMult.html">MatMult</a></span></code>(B,x,y) will compute the matrix-vector product of A transpose times x.</p>
</section>
<section id="other-matrix-operations">
<span id="sec-othermat"></span><h2>Other Matrix Operations<a class="headerlink" href="#other-matrix-operations" title="Link to this heading">#</a></h2>
<p>In many iterative calculations (for instance, in a nonlinear equations
solver), it is important for efficiency purposes to reuse the nonzero
structure of a matrix, rather than determining it anew every time the
matrix is generated. To retain a given matrix but reinitialize its
contents, one can employ</p>
<div class="highlight-c notranslate"><div class="highlight"><pre><span></span><span class="n"><a href="../manualpages/Mat/MatZeroEntries.html">MatZeroEntries</a></span><span class="p">(</span><span class="n"><a href="../manualpages/Mat/Mat.html">Mat</a></span><span class="w"> </span><span class="n">A</span><span class="p">);</span>
</pre></div>
</div>
<p>This routine will zero the matrix entries in the data structure but keep
all the data that indicates where the nonzeros are located. In this way
a new matrix assembly will be much less expensive, since no memory
allocations or copies will be needed. Of course, one can also explicitly
set selected matrix elements to zero by calling <code class="docutils notranslate"><span class="pre"><a href="../manualpages/Mat/MatSetValues.html">MatSetValues</a>()</span></code>.</p>
<p>By default, if new entries are made in locations where no nonzeros
previously existed, space will be allocated for the new entries. To
prevent the allocation of additional memory and simply discard those new
entries, one can use the option</p>
<div class="highlight-c notranslate"><div class="highlight"><pre><span></span><span class="n"><a href="../manualpages/Mat/MatSetOption.html">MatSetOption</a></span><span class="p">(</span><span class="n"><a href="../manualpages/Mat/Mat.html">Mat</a></span><span class="w"> </span><span class="n">A</span><span class="p">,</span><span class="n"><a href="../manualpages/Mat/MatOption.html">MAT_NEW_NONZERO_LOCATIONS</a></span><span class="p">,</span><span class="n"><a href="../manualpages/Sys/PETSC_FALSE.html">PETSC_FALSE</a></span><span class="p">);</span>
</pre></div>
</div>
<p>Once the matrix has been assembled, one can factor it numerically
without repeating the ordering or the symbolic factorization. This
option can save some computational time, although it does require that
the factorization is not done in-place.</p>
<p>In the numerical solution of elliptic partial differential equations, it
can be cumbersome to deal with Dirichlet boundary conditions. In
particular, one would like to assemble the matrix without regard to
boundary conditions and then at the end apply the Dirichlet boundary
conditions. In numerical analysis classes this process is usually
presented as moving the known boundary conditions to the right-hand side
and then solving a smaller linear system for the interior unknowns.
Unfortunately, implementing this requires extracting a large submatrix
from the original matrix and creating its corresponding data structures.
This process can be expensive in terms of both time and memory.</p>
<p>One simple way to deal with this difficulty is to replace those rows in
the matrix associated with known boundary conditions, by rows of the
identity matrix (or some scaling of it). This action can be done with
the command</p>
<div class="highlight-c notranslate"><div class="highlight"><pre><span></span><span class="n"><a href="../manualpages/Mat/MatZeroRows.html">MatZeroRows</a></span><span class="p">(</span><span class="n"><a href="../manualpages/Mat/Mat.html">Mat</a></span><span class="w"> </span><span class="n">A</span><span class="p">,</span><span class="n"><a href="../manualpages/Sys/PetscInt.html">PetscInt</a></span><span class="w"> </span><span class="n">numRows</span><span class="p">,</span><span class="n"><a href="../manualpages/Sys/PetscInt.html">PetscInt</a></span><span class="w"> </span><span class="n">rows</span><span class="p">[],</span><span class="n"><a href="../manualpages/Sys/PetscScalar.html">PetscScalar</a></span><span class="w"> </span><span class="n">diag_value</span><span class="p">,</span><span class="n"><a href="../manualpages/Vec/Vec.html">Vec</a></span><span class="w"> </span><span class="n">x</span><span class="p">,</span><span class="n"><a href="../manualpages/Vec/Vec.html">Vec</a></span><span class="w"> </span><span class="n">b</span><span class="p">),</span>
</pre></div>
</div>
<p>or equivalently,</p>
<div class="highlight-c notranslate"><div class="highlight"><pre><span></span><span class="n"><a href="../manualpages/Mat/MatZeroRowsIS.html">MatZeroRowsIS</a></span><span class="p">(</span><span class="n"><a href="../manualpages/Mat/Mat.html">Mat</a></span><span class="w"> </span><span class="n">A</span><span class="p">,</span><span class="n"><a href="../manualpages/IS/IS.html">IS</a></span><span class="w"> </span><span class="n">rows</span><span class="p">,</span><span class="n"><a href="../manualpages/Sys/PetscScalar.html">PetscScalar</a></span><span class="w"> </span><span class="n">diag_value</span><span class="p">,</span><span class="n"><a href="../manualpages/Vec/Vec.html">Vec</a></span><span class="w"> </span><span class="n">x</span><span class="p">,</span><span class="n"><a href="../manualpages/Vec/Vec.html">Vec</a></span><span class="w"> </span><span class="n">b</span><span class="p">);</span>
</pre></div>
</div>
<p>For sparse matrices this removes the data structures for certain rows of
the matrix. If the pointer <code class="docutils notranslate"><span class="pre">diag_value</span></code> is <code class="docutils notranslate"><span class="pre">NULL</span></code>, it even removes
the diagonal entry. If the pointer is not null, it uses that given value
at the pointer location in the diagonal entry of the eliminated rows.</p>
<p>One nice feature of this approach is that when solving a nonlinear
problem such that at each iteration the Dirichlet boundary conditions
are in the same positions and the matrix retains the same nonzero
structure, the user can call <code class="docutils notranslate"><span class="pre"><a href="../manualpages/Mat/MatZeroRows.html">MatZeroRows</a>()</span></code> in the first iteration.
Then, before generating the matrix in the second iteration the user
should call</p>
<div class="highlight-c notranslate"><div class="highlight"><pre><span></span><span class="n"><a href="../manualpages/Mat/MatSetOption.html">MatSetOption</a></span><span class="p">(</span><span class="n"><a href="../manualpages/Mat/Mat.html">Mat</a></span><span class="w"> </span><span class="n">A</span><span class="p">,</span><span class="n"><a href="../manualpages/Mat/MatOption.html">MAT_NEW_NONZERO_LOCATIONS</a></span><span class="p">,</span><span class="n"><a href="../manualpages/Sys/PETSC_FALSE.html">PETSC_FALSE</a></span><span class="p">);</span>
</pre></div>
</div>
<p>From that point, no new values will be inserted into those (boundary)
rows of the matrix.</p>
<p>The functions <code class="docutils notranslate"><span class="pre"><a href="../manualpages/Mat/MatZeroRowsLocal.html">MatZeroRowsLocal</a>()</span></code> and <code class="docutils notranslate"><span class="pre"><a href="../manualpages/Mat/MatZeroRowsLocalIS.html">MatZeroRowsLocalIS</a>()</span></code> can
also be used if for each process one provides the Dirichlet locations in
the local numbering of the matrix. A drawback of <code class="docutils notranslate"><span class="pre"><a href="../manualpages/Mat/MatZeroRows.html">MatZeroRows</a>()</span></code> is
that it destroys the symmetry of a matrix. Thus one can use</p>
<div class="highlight-c notranslate"><div class="highlight"><pre><span></span><span class="n"><a href="../manualpages/Mat/MatZeroRowsColumns.html">MatZeroRowsColumns</a></span><span class="p">(</span><span class="n"><a href="../manualpages/Mat/Mat.html">Mat</a></span><span class="w"> </span><span class="n">A</span><span class="p">,</span><span class="n"><a href="../manualpages/Sys/PetscInt.html">PetscInt</a></span><span class="w"> </span><span class="n">numRows</span><span class="p">,</span><span class="n"><a href="../manualpages/Sys/PetscInt.html">PetscInt</a></span><span class="w"> </span><span class="n">rows</span><span class="p">[],</span><span class="n"><a href="../manualpages/Sys/PetscScalar.html">PetscScalar</a></span><span class="w"> </span><span class="n">diag_value</span><span class="p">,</span><span class="n"><a href="../manualpages/Vec/Vec.html">Vec</a></span><span class="w"> </span><span class="n">x</span><span class="p">,</span><span class="n"><a href="../manualpages/Vec/Vec.html">Vec</a></span><span class="w"> </span><span class="n">b</span><span class="p">),</span>
</pre></div>
</div>
<p>or equivalently,</p>
<div class="highlight-c notranslate"><div class="highlight"><pre><span></span><span class="n"><a href="../manualpages/Mat/MatZeroRowsColumnsIS.html">MatZeroRowsColumnsIS</a></span><span class="p">(</span><span class="n"><a href="../manualpages/Mat/Mat.html">Mat</a></span><span class="w"> </span><span class="n">A</span><span class="p">,</span><span class="n"><a href="../manualpages/IS/IS.html">IS</a></span><span class="w"> </span><span class="n">rows</span><span class="p">,</span><span class="n"><a href="../manualpages/Sys/PetscScalar.html">PetscScalar</a></span><span class="w"> </span><span class="n">diag_value</span><span class="p">,</span><span class="n"><a href="../manualpages/Vec/Vec.html">Vec</a></span><span class="w"> </span><span class="n">x</span><span class="p">,</span><span class="n"><a href="../manualpages/Vec/Vec.html">Vec</a></span><span class="w"> </span><span class="n">b</span><span class="p">);</span>
</pre></div>
</div>
<p>Note that with all of these for a given assembled matrix it can be only
called once to update the x and b vector. It cannot be used if one
wishes to solve multiple right-hand side problems for the same matrix
since the matrix entries needed for updating the b vector are removed in
its first use.</p>
<p>Once the zeroed rows are removed the new matrix has possibly many rows
with only a diagonal entry affecting the parallel load balancing. The
<code class="docutils notranslate"><span class="pre"><a href="../manualpages/PC/PCREDISTRIBUTE.html">PCREDISTRIBUTE</a></span></code> preconditioner removes all the zeroed rows (and
associated columns and adjusts the right-hand side based on the removed
columns) and then rebalances the resulting rows of smaller matrix across
the processes. Thus one can use <code class="docutils notranslate"><span class="pre"><a href="../manualpages/Mat/MatZeroRows.html">MatZeroRows</a>()</span></code> to set the Dirichlet
points and then solve with the preconditioner <code class="docutils notranslate"><span class="pre"><a href="../manualpages/PC/PCREDISTRIBUTE.html">PCREDISTRIBUTE</a></span></code>. Note
if the original matrix was symmetric the smaller solved matrix will also
be symmetric.</p>
<p>Another matrix routine of interest is</p>
<div class="highlight-c notranslate"><div class="highlight"><pre><span></span><span class="n"><a href="../manualpages/Mat/MatConvert.html">MatConvert</a></span><span class="p">(</span><span class="n"><a href="../manualpages/Mat/Mat.html">Mat</a></span><span class="w"> </span><span class="n">mat</span><span class="p">,</span><span class="n"><a href="../manualpages/Mat/MatType.html">MatType</a></span><span class="w"> </span><span class="n">newtype</span><span class="p">,</span><span class="n"><a href="../manualpages/Mat/Mat.html">Mat</a></span><span class="w"> </span><span class="o">*</span><span class="n">M</span><span class="p">)</span>
</pre></div>
</div>
<p>which converts the matrix <code class="docutils notranslate"><span class="pre">mat</span></code> to new matrix, <code class="docutils notranslate"><span class="pre">M</span></code>, that has either
the same or different format. Set <code class="docutils notranslate"><span class="pre">newtype</span></code> to <code class="docutils notranslate"><span class="pre">MATSAME</span></code> to copy the
matrix, keeping the same matrix format. See
<code class="docutils notranslate"><span class="pre">$PETSC_DIR/include/petscmat.h</span></code> (<a href="../include/petscmat.h.html">source</a>)
for other available matrix types; standard ones are <code class="docutils notranslate"><span class="pre"><a href="../manualpages/Mat/MATSEQDENSE.html">MATSEQDENSE</a></span></code>,
<code class="docutils notranslate"><span class="pre"><a href="../manualpages/Mat/MATSEQAIJ.html">MATSEQAIJ</a></span></code>, <code class="docutils notranslate"><span class="pre"><a href="../manualpages/Mat/MATMPIAIJ.html">MATMPIAIJ</a></span></code>, <code class="docutils notranslate"><span class="pre"><a href="../manualpages/Mat/MATSEQBAIJ.html">MATSEQBAIJ</a></span></code> and <code class="docutils notranslate"><span class="pre"><a href="../manualpages/Mat/MATMPIBAIJ.html">MATMPIBAIJ</a></span></code>.</p>
<p>In certain applications it may be necessary for application codes to
directly access elements of a matrix. This may be done by using the the
command (for local rows only)</p>
<div class="highlight-c notranslate"><div class="highlight"><pre><span></span><span class="n"><a href="../manualpages/Mat/MatGetRow.html">MatGetRow</a></span><span class="p">(</span><span class="n"><a href="../manualpages/Mat/Mat.html">Mat</a></span><span class="w"> </span><span class="n">A</span><span class="p">,</span><span class="n"><a href="../manualpages/Sys/PetscInt.html">PetscInt</a></span><span class="w"> </span><span class="n">row</span><span class="p">,</span><span class="w"> </span><span class="n"><a href="../manualpages/Sys/PetscInt.html">PetscInt</a></span><span class="w"> </span><span class="o">*</span><span class="n">ncols</span><span class="p">,</span><span class="k">const</span><span class="w"> </span><span class="n"><a href="../manualpages/Sys/PetscInt.html">PetscInt</a></span><span class="w"> </span><span class="p">(</span><span class="o">*</span><span class="n">cols</span><span class="p">)[],</span><span class="k">const</span><span class="w"> </span><span class="n"><a href="../manualpages/Sys/PetscScalar.html">PetscScalar</a></span><span class="w"> </span><span class="p">(</span><span class="o">*</span><span class="n">vals</span><span class="p">)[]);</span>
</pre></div>
</div>
<p>The argument <code class="docutils notranslate"><span class="pre">ncols</span></code> returns the number of nonzeros in that row, while
<code class="docutils notranslate"><span class="pre">cols</span></code> and <code class="docutils notranslate"><span class="pre">vals</span></code> returns the column indices (with indices starting
at zero) and values in the row. If only the column indices are needed
(and not the corresponding matrix elements), one can use <code class="docutils notranslate"><span class="pre">NULL</span></code> for
the <code class="docutils notranslate"><span class="pre">vals</span></code> argument. Similarly, one can use <code class="docutils notranslate"><span class="pre">NULL</span></code> for the <code class="docutils notranslate"><span class="pre">cols</span></code>
argument. The user can only examine the values extracted with
<code class="docutils notranslate"><span class="pre"><a href="../manualpages/Mat/MatGetRow.html">MatGetRow</a>()</span></code>; the values <em>cannot</em> be altered. To change the matrix
entries, one must use <code class="docutils notranslate"><span class="pre"><a href="../manualpages/Mat/MatSetValues.html">MatSetValues</a>()</span></code>.</p>
<p>Once the user has finished using a row, he or she <em>must</em> call</p>
<div class="highlight-c notranslate"><div class="highlight"><pre><span></span><span class="n"><a href="../manualpages/Mat/MatRestoreRow.html">MatRestoreRow</a></span><span class="p">(</span><span class="n"><a href="../manualpages/Mat/Mat.html">Mat</a></span><span class="w"> </span><span class="n">A</span><span class="p">,</span><span class="n"><a href="../manualpages/Sys/PetscInt.html">PetscInt</a></span><span class="w"> </span><span class="n">row</span><span class="p">,</span><span class="n"><a href="../manualpages/Sys/PetscInt.html">PetscInt</a></span><span class="w"> </span><span class="o">*</span><span class="n">ncols</span><span class="p">,</span><span class="n"><a href="../manualpages/Sys/PetscInt.html">PetscInt</a></span><span class="w"> </span><span class="o">**</span><span class="n">cols</span><span class="p">,</span><span class="n"><a href="../manualpages/Sys/PetscScalar.html">PetscScalar</a></span><span class="w"> </span><span class="o">**</span><span class="n">vals</span><span class="p">);</span>
</pre></div>
</div>
<p>to free any space that was allocated during the call to <code class="docutils notranslate"><span class="pre"><a href="../manualpages/Mat/MatGetRow.html">MatGetRow</a>()</span></code>.</p>
</section>
<section id="symbolic-and-numeric-stages-in-sparse-matrix-operations">
<span id="sec-symbolic-numeric"></span><h2>Symbolic and Numeric Stages in Sparse Matrix Operations<a class="headerlink" href="#symbolic-and-numeric-stages-in-sparse-matrix-operations" title="Link to this heading">#</a></h2>
<p>Many sparse matrix operations can be optimized by dividing the computation into two stages: a symbolic stage that
creates any required data structures and does all the computations that do not require the matrices’ numerical values followed by one or more uses of a
numerical stage that use the symbolically computed information. Examples of such operations include <code class="docutils notranslate"><span class="pre"><a href="../manualpages/Mat/MatTranspose.html">MatTranspose</a>()</span></code>, <code class="docutils notranslate"><span class="pre"><a href="../manualpages/Mat/MatCreateSubMatrices.html">MatCreateSubMatrices</a>()</span></code>,
<code class="docutils notranslate"><span class="pre"><a href="../manualpages/Mat/MatCholeskyFactorSymbolic.html">MatCholeskyFactorSymbolic</a>()</span></code>, and <code class="docutils notranslate"><span class="pre"><a href="../manualpages/Mat/MatCholeskyFactorNumeric.html">MatCholeskyFactorNumeric</a>()</span></code>.
PETSc uses two different API’s to take advantage of these optimizations.</p>
<p>The first approach explicitly divides the computation in the API. This approach is used, for example, with <code class="docutils notranslate"><span class="pre"><a href="../manualpages/Mat/MatCholeskyFactorSymbolic.html">MatCholeskyFactorSymbolic</a>()</span></code>, <code class="docutils notranslate"><span class="pre"><a href="../manualpages/Mat/MatCholeskyFactorNumeric.html">MatCholeskyFactorNumeric</a>()</span></code>.
The caller can take advantage of their knowledge of changes in the nonzero structure of the sparse matrices to call the appropriate routines as needed. In fact, they can
use <code class="docutils notranslate"><span class="pre"><a href="../manualpages/Mat/MatGetNonzeroState.html">MatGetNonzeroState</a>()</span></code> to determine if a new symbolic computation is needed. The drawback of this approach is that the caller of these routines has to
manage the creation of new matrices when the nonzero structure changes.</p>
<p>The second approach, as exemplified by <code class="docutils notranslate"><span class="pre"><a href="../manualpages/Mat/MatTranspose.html">MatTranspose</a>()</span></code>, does not expose the two stages explicit in the API, instead a flag, <code class="docutils notranslate"><span class="pre"><a href="../manualpages/Mat/MatReuse.html">MatReuse</a></span></code> is passed through the
API to indicate if a symbolic data structure is already available or needs to be computed. Thus <code class="docutils notranslate"><span class="pre"><a href="../manualpages/Mat/MatTranspose.html">MatTranspose</a>(A,<a href="../manualpages/Mat/MatReuse.html">MAT_INITIAL_MATRIX</a>,&amp;B)</span></code> is called first, then
<code class="docutils notranslate"><span class="pre"><a href="../manualpages/Mat/MatTranspose.html">MatTranspose</a>(A,<a href="../manualpages/Mat/MatReuse.html">MAT_REUSE_MATRIX</a>,&amp;B)</span></code> can be called repeatedly with new numerical values in the A matrix. In theory, if the nonzero structure of A changes, the
symbolic computations for B could be redone automatically inside the same B matrix when there is a change in the nonzero state of the A matrix. In practice, in PETSc, the
<code class="docutils notranslate"><span class="pre"><a href="../manualpages/Mat/MatReuse.html">MAT_REUSE_MATRIX</a></span></code> for most PETSc routines only works if the nonzero structure does not change and the code may crash otherwise. The advantage of this approach
(when the nonzero structure changes are handled correctly) is that the calling code does not need to keep track of the nonzero state of the matrices; everything
“just works”. However, the caller must still know when it is the first call to the routine so the flag <code class="docutils notranslate"><span class="pre"><a href="../manualpages/Mat/MatReuse.html">MAT_INITIAL_MATRIX</a></span></code> is being used. If the underlying implementation language supported detecting a yet to be initialized variable at run time, the <code class="docutils notranslate"><span class="pre"><a href="../manualpages/Mat/MatReuse.html">MatReuse</a></span></code> flag would not be need.</p>
<p>PETSc uses two approaches because the same programming problem was solved with two different ways during PETSc’s early development.
A better model would combine both approaches; an explicit
separation of the stages and a unified operation that internally utilized the two stages appropriately and also handled changes to the nonzero structure. Code could be simplified in many places with this approach, in most places the use of the unified API would replace the use of the separate stages.</p>
<p>See <a class="reference internal" href="advanced.html#sec-matsub"><span class="std std-ref">Extracting Submatrices</span></a> and <a class="reference internal" href="advanced.html#sec-matmatproduct"><span class="std std-ref">Matrix-Matrix Products</span></a>.</p>
</section>
<section id="graph-operations">
<span id="sec-graph"></span><h2>Graph Operations<a class="headerlink" href="#graph-operations" title="Link to this heading">#</a></h2>
<p>PETSc has four families of graph operations that treat sparse <code class="docutils notranslate"><span class="pre"><a href="../manualpages/Mat/Mat.html">Mat</a></span></code> as representing graphs.</p>
<table class="table table-center">
<thead>
<tr class="row-odd"><th class="head"><p>Operation</p></th>
<th class="head"><p>Type</p></th>
<th class="head"><p>Available methods</p></th>
<th class="head"><p>User guide</p></th>
</tr>
</thead>
<tbody>
<tr class="row-even"><td><p>Ordering to reduce fill</p></td>
<td><p>N/A</p></td>
<td><p><code class="docutils notranslate"><span class="pre"><a href="../manualpages/Mat/MatOrderingType.html">MatOrderingType</a></span></code></p></td>
<td><p><a class="reference internal" href="advanced.html#sec-matfactor"><span class="std std-ref">Matrix Factorization</span></a></p></td>
</tr>
<tr class="row-odd"><td><p>Partitioning for parallelism</p></td>
<td><p><code class="docutils notranslate"><span class="pre"><a href="../manualpages/Mat/MatPartitioning.html">MatPartitioning</a></span></code></p></td>
<td><p><code class="docutils notranslate"><span class="pre"><a href="../manualpages/Mat/MatPartitioningType.html">MatPartitioningType</a></span></code></p></td>
<td><p><a class="reference internal" href="#sec-partitioning"><span class="std std-ref">Partitioning</span></a></p></td>
</tr>
<tr class="row-even"><td><p>Coloring for parallelism or Jacobians</p></td>
<td><p><code class="docutils notranslate"><span class="pre"><a href="../manualpages/Mat/MatColoring.html">MatColoring</a></span></code></p></td>
<td><p><code class="docutils notranslate"><span class="pre"><a href="../manualpages/Mat/MatColoringType.html">MatColoringType</a></span></code></p></td>
<td><p><a class="reference internal" href="snes.html#sec-fdmatrix"><span class="std std-ref">Finite Difference Jacobian Approximations</span></a></p></td>
</tr>
<tr class="row-odd"><td><p>Coarsening for multigrid</p></td>
<td><p><code class="docutils notranslate"><span class="pre"><a href="../manualpages/Mat/MatCoarsen.html">MatCoarsen</a></span></code></p></td>
<td><p><code class="docutils notranslate"><span class="pre"><a href="../manualpages/Mat/MatCoarsenType.html">MatCoarsenType</a></span></code></p></td>
<td><p><a class="reference internal" href="ksp.html#sec-amg"><span class="std std-ref">Algebraic Multigrid (AMG) Preconditioners</span></a></p></td>
</tr>
</tbody>
</table>
</section>
<section id="partitioning">
<span id="sec-partitioning"></span><h2>Partitioning<a class="headerlink" href="#partitioning" title="Link to this heading">#</a></h2>
<p>For almost all unstructured grid computation, the distribution of
portions of the grid across the process’s work load and memory can have
a very large impact on performance. In most PDE calculations the grid
partitioning and distribution across the processes can (and should) be
done in a “pre-processing” step before the numerical computations.
However, this does not mean it need be done in a separate, sequential
program; rather, it should be done before one sets up the parallel grid
data structures in the actual program. PETSc provides an interface to
the ParMETIS (developed by George Karypis; see
<a class="reference external" href="https://petsc.org/release/install/">the PETSc installation instructions</a>
for directions on installing PETSc to use ParMETIS) to allow the
partitioning to be done in parallel. PETSc does not currently provide
directly support for dynamic repartitioning, load balancing by migrating
matrix entries between processes, etc. For problems that require mesh
refinement, PETSc uses the “rebuild the data structure” approach, as
opposed to the “maintain dynamic data structures that support the
insertion/deletion of additional vector and matrix rows and columns
entries” approach.</p>
<p>Partitioning in PETSc is organized around the <code class="docutils notranslate"><span class="pre"><a href="../manualpages/Mat/MatPartitioning.html">MatPartitioning</a></span></code>
object. One first creates a parallel matrix that contains the
connectivity information about the grid (or other graph-type object)
that is to be partitioned. This is done with the command</p>
<div class="highlight-c notranslate"><div class="highlight"><pre><span></span><span class="n"><a href="../manualpages/Mat/MatCreateMPIAdj.html">MatCreateMPIAdj</a></span><span class="p">(</span><span class="n"><a href="../manualpages/Sys/MPI_Comm.html">MPI_Comm</a></span><span class="w"> </span><span class="n">comm</span><span class="p">,</span><span class="kt">int</span><span class="w"> </span><span class="n">mlocal</span><span class="p">,</span><span class="n"><a href="../manualpages/Sys/PetscInt.html">PetscInt</a></span><span class="w"> </span><span class="n">n</span><span class="p">,</span><span class="k">const</span><span class="w"> </span><span class="n"><a href="../manualpages/Sys/PetscInt.html">PetscInt</a></span><span class="w"> </span><span class="n">ia</span><span class="p">[],</span><span class="k">const</span><span class="w"> </span><span class="n"><a href="../manualpages/Sys/PetscInt.html">PetscInt</a></span><span class="w"> </span><span class="n">ja</span><span class="p">[],</span><span class="n"><a href="../manualpages/Sys/PetscInt.html">PetscInt</a></span><span class="w"> </span><span class="o">*</span><span class="n">weights</span><span class="p">,</span><span class="n"><a href="../manualpages/Mat/Mat.html">Mat</a></span><span class="w"> </span><span class="o">*</span><span class="n">Adj</span><span class="p">);</span>
</pre></div>
</div>
<p>The argument <code class="docutils notranslate"><span class="pre">mlocal</span></code> indicates the number of rows of the graph being
provided by the given process, <code class="docutils notranslate"><span class="pre">n</span></code> is the total number of columns;
equal to the sum of all the <code class="docutils notranslate"><span class="pre">mlocal</span></code>. The arguments <code class="docutils notranslate"><span class="pre">ia</span></code> and <code class="docutils notranslate"><span class="pre">ja</span></code>
are the row pointers and column pointers for the given rows; these are
the usual format for parallel compressed sparse row storage, using
indices starting at 0, <em>not</em> 1.</p>
<figure class="align-default" id="fig-usg">
<img alt="Numbering on Simple Unstructured Grid" src="../_images/usg.png" />
<figcaption>
<p><span class="caption-number">Fig. 6 </span><span class="caption-text">Numbering on Simple Unstructured Grid</span><a class="headerlink" href="#fig-usg" title="Link to this image">#</a></p>
</figcaption>
</figure>
<p>This, of course, assumes that one has already distributed the grid
(graph) information among the processes. The details of this initial
distribution is not important; it could be simply determined by
assigning to the first process the first <span class="math">\(n_0\)</span> nodes from a file,
the second process the next <span class="math">\(n_1\)</span> nodes, etc.</p>
<p>For example, we demonstrate the form of the <code class="docutils notranslate"><span class="pre">ia</span></code> and <code class="docutils notranslate"><span class="pre">ja</span></code> for a
triangular grid where we</p>
<ol class="arabic simple">
<li><p>partition by element (triangle)</p></li>
</ol>
<ul class="simple">
<li><p>Process 0: <code class="docutils notranslate"><span class="pre">mlocal</span> <span class="pre">=</span> <span class="pre">2</span></code>, <code class="docutils notranslate"><span class="pre">n</span> <span class="pre">=</span> <span class="pre">4</span></code>, <code class="docutils notranslate"><span class="pre">ja</span> <span class="pre">=``{2,3,</span> <span class="pre">3}</span></code>,
<code class="docutils notranslate"><span class="pre">ia</span> <span class="pre">=</span></code> <code class="docutils notranslate"><span class="pre">{0,2,3}</span></code></p></li>
<li><p>Process 1: <code class="docutils notranslate"><span class="pre">mlocal</span> <span class="pre">=</span> <span class="pre">2</span></code>, <code class="docutils notranslate"><span class="pre">n</span> <span class="pre">=</span> <span class="pre">4</span></code>, <code class="docutils notranslate"><span class="pre">ja</span> <span class="pre">=``{0,</span> <span class="pre">0,1}</span></code>,
<code class="docutils notranslate"><span class="pre">ia</span> <span class="pre">=``{0,1,3}</span></code></p></li>
</ul>
<p>Note that elements are not connected to themselves and we only indicate
edge connections (in some contexts single vertex connections between
elements may also be included). We use a space above to denote the
transition between rows in the matrix; and</p>
<ol class="arabic simple" start="2">
<li><p>partition by vertex.</p></li>
</ol>
<ul class="simple">
<li><p>Process 0: <code class="docutils notranslate"><span class="pre">mlocal</span> <span class="pre">=</span> <span class="pre">3</span></code>, <code class="docutils notranslate"><span class="pre">n</span> <span class="pre">=</span> <span class="pre">6</span></code>,
<code class="docutils notranslate"><span class="pre">ja</span> <span class="pre">=``{3,4,</span> <span class="pre">4,5,</span> <span class="pre">3,4,5}</span></code>, <code class="docutils notranslate"><span class="pre">ia</span> <span class="pre">=``{0,</span> <span class="pre">2,</span> <span class="pre">4,</span> <span class="pre">7}</span></code></p></li>
<li><p>Process 1: <code class="docutils notranslate"><span class="pre">mlocal</span> <span class="pre">=</span> <span class="pre">3</span></code>, <code class="docutils notranslate"><span class="pre">n</span> <span class="pre">=</span> <span class="pre">6</span></code>,
<code class="docutils notranslate"><span class="pre">ja</span> <span class="pre">=``{0,2,</span> <span class="pre">4,</span> <span class="pre">0,1,2,3,5,</span> <span class="pre">1,2,4}</span></code>,
<code class="docutils notranslate"><span class="pre">ia</span> <span class="pre">=``{0,</span> <span class="pre">3,</span> <span class="pre">8,</span> <span class="pre">11}</span></code>.</p></li>
</ul>
<p>Once the connectivity matrix has been created the following code will
generate the renumbering required for the new partition</p>
<div class="highlight-c notranslate"><div class="highlight"><pre><span></span><span class="n"><a href="../manualpages/MatGraphOperations/MatPartitioningCreate.html">MatPartitioningCreate</a></span><span class="p">(</span><span class="n"><a href="../manualpages/Sys/MPI_Comm.html">MPI_Comm</a></span><span class="w"> </span><span class="n">comm</span><span class="p">,</span><span class="n"><a href="../manualpages/Mat/MatPartitioning.html">MatPartitioning</a></span><span class="w"> </span><span class="o">*</span><span class="n">part</span><span class="p">);</span>
<span class="n"><a href="../manualpages/MatGraphOperations/MatPartitioningSetAdjacency.html">MatPartitioningSetAdjacency</a></span><span class="p">(</span><span class="n"><a href="../manualpages/Mat/MatPartitioning.html">MatPartitioning</a></span><span class="w"> </span><span class="n">part</span><span class="p">,</span><span class="n"><a href="../manualpages/Mat/Mat.html">Mat</a></span><span class="w"> </span><span class="n">Adj</span><span class="p">);</span>
<span class="n"><a href="../manualpages/MatGraphOperations/MatPartitioningSetFromOptions.html">MatPartitioningSetFromOptions</a></span><span class="p">(</span><span class="n"><a href="../manualpages/Mat/MatPartitioning.html">MatPartitioning</a></span><span class="w"> </span><span class="n">part</span><span class="p">);</span>
<span class="n"><a href="../manualpages/MatGraphOperations/MatPartitioningApply.html">MatPartitioningApply</a></span><span class="p">(</span><span class="n"><a href="../manualpages/Mat/MatPartitioning.html">MatPartitioning</a></span><span class="w"> </span><span class="n">part</span><span class="p">,</span><span class="n"><a href="../manualpages/IS/IS.html">IS</a></span><span class="w"> </span><span class="o">*</span><span class="n">is</span><span class="p">);</span>
<span class="n"><a href="../manualpages/MatGraphOperations/MatPartitioningDestroy.html">MatPartitioningDestroy</a></span><span class="p">(</span><span class="n"><a href="../manualpages/Mat/MatPartitioning.html">MatPartitioning</a></span><span class="w"> </span><span class="o">*</span><span class="n">part</span><span class="p">);</span>
<span class="n"><a href="../manualpages/Mat/MatDestroy.html">MatDestroy</a></span><span class="p">(</span><span class="n"><a href="../manualpages/Mat/Mat.html">Mat</a></span><span class="w"> </span><span class="o">*</span><span class="n">Adj</span><span class="p">);</span>
<span class="n"><a href="../manualpages/IS/ISPartitioningToNumbering.html">ISPartitioningToNumbering</a></span><span class="p">(</span><span class="n"><a href="../manualpages/IS/IS.html">IS</a></span><span class="w"> </span><span class="n">is</span><span class="p">,</span><span class="n"><a href="../manualpages/IS/IS.html">IS</a></span><span class="w"> </span><span class="o">*</span><span class="n">isg</span><span class="p">);</span>
</pre></div>
</div>
<p>The resulting <code class="docutils notranslate"><span class="pre">isg</span></code> contains for each local node the new global number
of that node. The resulting <code class="docutils notranslate"><span class="pre">is</span></code> contains the new process number that
each local node has been assigned to.</p>
<p>Now that a new numbering of the nodes has been determined, one must
renumber all the nodes and migrate the grid information to the correct
process. The command</p>
<div class="highlight-c notranslate"><div class="highlight"><pre><span></span><span class="n"><a href="../manualpages/AO/AOCreateBasicIS.html">AOCreateBasicIS</a></span><span class="p">(</span><span class="n">isg</span><span class="p">,</span><span class="nb">NULL</span><span class="p">,</span><span class="o">&amp;</span><span class="n">ao</span><span class="p">);</span>
</pre></div>
</div>
<p>generates, see <a class="reference internal" href="vec.html#sec-ao"><span class="std std-ref">Application Orderings</span></a>, an <code class="docutils notranslate"><span class="pre"><a href="../manualpages/AO/AO.html">AO</a></span></code> object that can be
used in conjunction with the <code class="docutils notranslate"><span class="pre">is</span></code> and <code class="docutils notranslate"><span class="pre">isg</span></code> to move the relevant
grid information to the correct process and renumber the nodes etc. In
this context, the new ordering is the “application” ordering so
<code class="docutils notranslate"><span class="pre"><a href="../manualpages/AO/AOPetscToApplication.html">AOPetscToApplication</a>()</span></code> converts old global indices to new global
indices and <code class="docutils notranslate"><span class="pre"><a href="../manualpages/AO/AOApplicationToPetsc.html">AOApplicationToPetsc</a>()</span></code> converts new global indices back
to old global indices.</p>
<p>PETSc does not currently provide tools that completely manage the
migration and node renumbering, since it will be dependent on the
particular data structure you use to store the grid information and the
type of grid information that you need for your application. We do plan
to include more support for this in the future, but designing the
appropriate general user interface and providing a scalable
implementation that can be used for a wide variety of different grids
requires a great deal of time.</p>
<p>See <a class="reference internal" href="snes.html#sec-fdmatrix"><span class="std std-ref">Finite Difference Jacobian Approximations</span></a> and <a class="reference internal" href="advanced.html#sec-matfactor"><span class="std std-ref">Matrix Factorization</span></a> for discussions on performing graph coloring and computing graph reorderings to
reduce fill in sparse matrix factorizations.</p>
<div class="docutils container" id="id1">
<div role="list" class="citation-list">
<div class="citation" id="id2238" role="doc-biblioentry">
<span class="label"><span class="fn-bracket">[</span><a role="doc-backlink" href="#id8">ref-EM17a</a><span class="fn-bracket">]</span></span>
<p>Jennifer B Erway and Roummel F Marcia. On solving large-scale limited-memory quasi-newton equations. <em>Linear Algebra and its Applications</em>, 515:196–225, 2017.</p>
</div>
<div class="citation" id="id3945" role="doc-biblioentry">
<span class="label"><span class="fn-bracket">[</span>ref-EM17b<span class="fn-bracket">]</span></span>
<span class="backrefs">(<a role="doc-backlink" href="#id5">1</a>,<a role="doc-backlink" href="#id7">2</a>,<a role="doc-backlink" href="#id3">3</a>)</span>
<p>Jennifer B Erway and Roummel F Marcia. On solving large-scale limited-memory quasi-Newton equations. <em>Linear Algebra and its Applications</em>, 515:196–225, 2017.</p>
</div>
<div class="citation" id="id2905" role="doc-biblioentry">
<span class="label"><span class="fn-bracket">[</span>ref-GL89<span class="fn-bracket">]</span></span>
<span class="backrefs">(<a role="doc-backlink" href="#id1">1</a>,<a role="doc-backlink" href="#id2">2</a>)</span>
<p>J. C. Gilbert and C. Lemarechal. Some numerical experiments with variable-storage quasi-newton algorithms. <em>Mathematical Programming</em>, 45:407–434, 1989.</p>
</div>
<div class="citation" id="id2239" role="doc-biblioentry">
<span class="label"><span class="fn-bracket">[</span>ref-Gri12<span class="fn-bracket">]</span></span>
<span class="backrefs">(<a role="doc-backlink" href="#id1">1</a>,<a role="doc-backlink" href="#id2">2</a>)</span>
<p>Andreas Griewank. Broyden updating, the good and the bad! <em>Optimization Stories, Documenta Mathematica. Extra Volume: Optimization Stories</em>, pages 301–315, 2012.</p>
</div>
<div class="citation" id="id3887" role="doc-biblioentry">
<span class="label"><span class="fn-bracket">[</span>ref-NW06<span class="fn-bracket">]</span></span>
<span class="backrefs">(<a role="doc-backlink" href="#id3">1</a>,<a role="doc-backlink" href="#id4">2</a>,<a role="doc-backlink" href="#id6">3</a>)</span>
<p>Jorge Nocedal and Stephen Wright. <em>Numerical optimization</em>. Springer Science &amp; Business Media, 2006.</p>
</div>
</div>
</div>
</section>
</section>


                </article>
              
              
              
              
              
                <footer class="prev-next-footer">
                  
<div class="prev-next-area">
    <a class="left-prev"
       href="vec.html"
       title="previous page">
      <i class="fa-solid fa-angle-left"></i>
      <div class="prev-next-info">
        <p class="prev-next-subtitle">previous</p>
        <p class="prev-next-title">Vectors and Parallel Data</p>
      </div>
    </a>
    <a class="right-next"
       href="ksp.html"
       title="next page">
      <div class="prev-next-info">
        <p class="prev-next-subtitle">next</p>
        <p class="prev-next-title">KSP: Linear System Solvers</p>
      </div>
      <i class="fa-solid fa-angle-right"></i>
    </a>
</div>
                </footer>
              
            </div>
            
            
              
                <div class="bd-sidebar-secondary bd-toc"><div class="sidebar-secondary-items sidebar-secondary__inner">


  <div class="sidebar-secondary-item">
<div
    id="pst-page-navigation-heading-2"
    class="page-toc tocsection onthispage">
    <i class="fa-solid fa-list"></i> On this page
  </div>
  <nav class="bd-toc-nav page-toc" aria-labelledby="pst-page-navigation-heading-2">
    <ul class="visible nav section-nav flex-column">
<li class="toc-h2 nav-item toc-entry"><a class="reference internal nav-link" href="#creating-matrices">Creating matrices</a></li>
<li class="toc-h2 nav-item toc-entry"><a class="reference internal nav-link" href="#low-level-matrix-creation-routines">Low-level matrix creation routines</a></li>
<li class="toc-h2 nav-item toc-entry"><a class="reference internal nav-link" href="#assembling-putting-values-into-matrices">Assembling (putting values into) matrices</a><ul class="nav section-nav flex-column">
<li class="toc-h3 nav-item toc-entry"><a class="reference internal nav-link" href="#matrix-and-vector-layouts-and-storage-locations">Matrix and Vector Layouts and Storage Locations</a></li>
<li class="toc-h3 nav-item toc-entry"><a class="reference internal nav-link" href="#sparse-matrices">Sparse Matrices</a><ul class="nav section-nav flex-column">
<li class="toc-h4 nav-item toc-entry"><a class="reference internal nav-link" href="#sequential-aij-sparse-matrices">Sequential AIJ Sparse Matrices</a></li>
<li class="toc-h4 nav-item toc-entry"><a class="reference internal nav-link" href="#preallocation-of-memory-for-sequential-aij-sparse-matrices">Preallocation of Memory for Sequential AIJ Sparse Matrices</a></li>
<li class="toc-h4 nav-item toc-entry"><a class="reference internal nav-link" href="#parallel-aij-sparse-matrices">Parallel AIJ Sparse Matrices</a></li>
<li class="toc-h4 nav-item toc-entry"><a class="reference internal nav-link" href="#preallocation-of-memory-for-parallel-aij-sparse-matrices">Preallocation of Memory for Parallel AIJ Sparse Matrices</a></li>
<li class="toc-h4 nav-item toc-entry"><a class="reference internal nav-link" href="#limited-memory-variable-metric-lmvm-matrices">Limited-Memory Variable Metric (LMVM) Matrices</a></li>
</ul>
</li>
<li class="toc-h3 nav-item toc-entry"><a class="reference internal nav-link" href="#dense-matrices">Dense Matrices</a></li>
<li class="toc-h3 nav-item toc-entry"><a class="reference internal nav-link" href="#block-matrices">Block Matrices</a></li>
</ul>
</li>
<li class="toc-h2 nav-item toc-entry"><a class="reference internal nav-link" href="#basic-matrix-operations">Basic Matrix Operations</a></li>
<li class="toc-h2 nav-item toc-entry"><a class="reference internal nav-link" href="#application-specific-custom-matrices">Application Specific Custom Matrices</a></li>
<li class="toc-h2 nav-item toc-entry"><a class="reference internal nav-link" href="#transposes-of-matrices">Transposes of Matrices</a></li>
<li class="toc-h2 nav-item toc-entry"><a class="reference internal nav-link" href="#other-matrix-operations">Other Matrix Operations</a></li>
<li class="toc-h2 nav-item toc-entry"><a class="reference internal nav-link" href="#symbolic-and-numeric-stages-in-sparse-matrix-operations">Symbolic and Numeric Stages in Sparse Matrix Operations</a></li>
<li class="toc-h2 nav-item toc-entry"><a class="reference internal nav-link" href="#graph-operations">Graph Operations</a></li>
<li class="toc-h2 nav-item toc-entry"><a class="reference internal nav-link" href="#partitioning">Partitioning</a></li>
</ul>
  </nav></div>

  <div class="sidebar-secondary-item">

  
  <div class="tocsection editthispage">
    <a href="https://gitlab.com/petsc/petsc/-/edit/release/doc/manual/mat.md">
      <i class="fa-solid fa-pencil"></i>
      
      
        
          Edit on GitLab
        
      
    </a>
  </div>
</div>

  <div class="sidebar-secondary-item">

  <div class="tocsection sourcelink">
    <a href="../_sources/manual/mat.md.txt">
      <i class="fa-solid fa-file-lines"></i> Show Source
    </a>
  </div>
</div>

</div></div>
              
            
          </div>
          <footer class="bd-footer-content">
            
          </footer>
        
      </main>
    </div>
  </div>
  
  <!-- Scripts loaded after <body> so the DOM is not blocked -->
  <script src="../_static/scripts/bootstrap.js?digest=bd9e20870c6007c4c509"></script>
<script src="../_static/scripts/pydata-sphinx-theme.js?digest=bd9e20870c6007c4c509"></script>

  <footer class="bd-footer">
<div class="bd-footer__inner bd-page-width">
  
    <div class="footer-items__start">
      
        <div class="footer-item">

  <p class="copyright">
    
      © Copyright 1991-2025, UChicago Argonne, LLC and the PETSc Development Team.
      <br/>
    
  </p>
</div>
      
        <div class="footer-item">

  <p class="sphinx-version">
    Created using <a href="https://www.sphinx-doc.org/">Sphinx</a> 7.3.7.
    <br/>
  </p>
</div>
      
    </div>
  
  
  
    <div class="footer-items__end">
      
        <div class="footer-item">
<p class="theme-version">
  Built with the <a href="https://pydata-sphinx-theme.readthedocs.io/en/stable/index.html">PyData Sphinx Theme</a> 0.15.1.
</p></div>
      
        <div class="footer-item"><p class="last-updated">
  Last updated on 2025-04-30T13:10:40-0500 (v3.23.1).
  <br/>
</p></div>
      
    </div>
  
</div>

  </footer>
  </body>
</html>