Squashed 'third_party/ceres/' changes from 399cda773..7f9cc571b

7f9cc571b Set CMAKE_CUDA_ARCHITECTURES depending on CUDAToolkit_VERSION
6a74af202 Remove a level of indirection when using CellInfo
f8c2994da Drop Ubuntu 20.04 and add Ubuntu 24.04 support
20954e693 Eliminate CUDA set but unused variable warning
57aba3ed0 Enable Apple linker library deduplication
1f15197be Drop macos-11 runner and support macos-14 instead
5de0fda0f Update Github actions
308a5bb43 Add missing include
715865101 Remove 32-bit MinGW from CI matrix
f71181a92 Remove remaining references to CXSparse
522210a08 Reuse macro to format version string
cd2dd06e9 Fix clang16 compiler warnings
1f2e6313a Add missing std qualifiers
125c06882 Link static cuda libs when ceres is build static
62c03d6ff Revert "Update Eigen::BDCSVD usage to comply with Eigen 3.4"
4027f6997 Update Eigen::BDCSVD usage to comply with Eigen 3.4
1b2ebabf5 Typo fixes in the documentation
2ffeb943a Fix typo in AutoDiffManifold comment and docs
da34da3da Remove CreateFakeBundleAdjustmentPartitionedJacobian
85b2c418a ClangTidy fixes
91773746b Simplify instantiation of cost functions and their functors
8b88a9ab4 Use C++17 Bessel functions
84436f791 Use native CMake TBB package configuration
556a56f21 Do not assume Homebrew usage
3fd2a72cc MinGW no longer provides a 32-bit METIS package
e776a6f1a Unbreak the Bazel build.
095d48392 Skip structure detection on preprocessor failure
a0876309a Fix typos in comments
774973731 Fix search on ceres-solver.org
85331393d Update docs for 2.2.0.
2120eae67 Optimize the computation of the LM diagonal in TinySolver
611b139b1 Fix Solver::Options::callbacks type in documentation
b652d3b4f Unbreak the windows build
b79c4d350 ClangTidy fixes
76af132d0 Update docs for 2.2.0RC3
dc7a85975 Single-threaded operations on small vectors
b379ab768 Remove MaxNumThreadsAvailable
354002f98 Schedule task in ParallerFor from the previous task
a9b3fcff4 Minor update to docs
5ccab18be Drop use of POSIX M_PI_2 and M_PI_4
4519b8d77 Drop use of POSIX M_PI
b83abdcb1 Add a default value for Solver::Summary::linear_solver_ordering_type
8d875a312 Fix checks for CUDA memory pools support
94335e3b9 More ClangTidy fixes
399395c4f Miscellaneous ClangTidy fixes
c8bed4b93 Update version_history for 2.2.0rc2
489339219 Rework MSVC warning suppression
0cea191d4 Move stream-ordered memory allocations
8e3b7d89e Fix a copy-pasta error
dc0bb8508 Various cleanups to the documentation
4588b0fbb Add an example for EvaluationCallback
0fc3fd4fc Add documentation for examples
e6b2f532b Parallelize PSE preconditioner
41672dff8 Add end-to-end BA tests for SCHUR_POWER_SERIES_EXPANSION
83ee376d8 Add an example for IterationCallback
371265094 Cleanup example code
59182a42c Update documentation
d4db6e6fe Fix typos in the documentation for EvaluationCallback
dffd8cd71 Add an accessor for the CostFunctor in DynamicAutoDiffCostFunction
bea247701 Add a missing include dir to the cuda kernels target.
18ea7d1c2 Runtime check for cudaMallocAsync support
a227045be Remove cuda-memcheck based tests
d10e786ca Remove an unused variable from CudaSparseMatrix
5a30cae58 Preparing for 2.2.0rc1
9cca67127 Enable compatibility with SuiteSparse 7.2.0
a1c02e8d3 Rework the Sphinx find module
a57e35bba Require at least CMake 3.16
863db948f Eliminate macOS sprintf warning
6f5342db6 Export ceres_cuda_kernels to project's build tree
d864d146f Add macOS 13 runner to Github workflows
01a23504d Add a workaround for CMake Policy CMP0148
de9cbde95 Work around MinGW32 manifold_test segfault
5e4b22f7f Update CudaSparseMatrix class
ed9921fc2 Fix Solver::Options in documentation
de62bf220 Two minor fixes
5f97455be Fix typos in documentation
5ba62abec Add CUDA support to windows CI builds
a98fdf582 Update CMakeLists.txt to fix Windows CUDA build
ec4907399 Fix block-sparse to crs conversion on windows
799ee91bb Fix check in CompressedRowJacobianWriter::CreateJacobian()
ee90f6cbf Detect large Jacobians and return failure instead of crashing.
310a252fb Deal with infinite initial cost correctly.
357482db7 Add NumTraits::max_digits10 for Jets
908b5b1b5 Fix type mismatch in documentation
75bacedf7 CUDA partitioned matrix view
fd6197ce0 Fixed GCC 13.1 compilation errors
df97a8f05 Improve support of older CUDA toolkit versions
085214ea7 Fix test CompressedRowSparseMatrix.Transpose
d880df09f Match new[] with delete[] in BSM
bdee4d617 Block-sparse to CRS conversion using block-structure
0f9de3daf Use page locked memory in BlockSparseMatrix
e7bd72d41 Permutation-based conversion from block-sparse to crs
abbc4e797 Explicitly compute number of non-zeros in row
96fdfd2e7 Implement tests for Euler conversion with jets
db1ebd3ff Work around lack of constexpr constructors for Jet
16a4fa04e Further Jet conversion fixes
92ad18b8a Fix a Jet conversion bug in rotation.h
a5e745d4e ClangTidy fixes
77ad8bb4e Change storage in BlockRandomAccessSparseMatrix
d340f81bd Clang Tidy fixes
54ad3dd03 Reorganize ParallelFor source files
ba360ab07 Change the value of BlockRandomAccessSparseMatrix::kMaxRowBlocks
0315c6ca9 Provide DynamicAutoDiffCostFunction deduction guide
3cdfae110 Replace Hertzberg mentions with citations
0af38a9fc Fix typos in documentation
4c969a6c1 Improve image of loss functions shape
b54f05b8e Add missing TukeyLoss to documentation
8bf4a2f42 Inexact check for ParallelAssign test
9cddce73a Explicit conversions from long to int in benchmarks (for num_threads)
5e787ab70 Using const ints in jacobian writers
f9bffbb6f Removing -Wshorten-64-to-32 warnings from examples (part 1)
e269b64f5 More ClangTidy fixes
f4eb768e0 Using int64 in file.cc. Fixing compilation error in array_utils.cc
74a0f0d24 Using int64_t for sizes and indexes in array utils
749a442d9 Clang-Tidy fixes
c4ba975ae Fix rotation_test.cc to work with older versions of Eigen
9602ed7b7 ClangFormat changes
79a554ffc Fix a bug in QuaternionRotatePoint.
f1113c08a Commenting unused parameters for better readibility
772d927e1 Replacing old style typedefs with new style usings
53df5ddcf Removing using std::...
e1ca3302a Increasing bazel timeout for heavy BA tests
f982d3071 Fixing bazel build
cb6b30662 Use hypot to compute the L^2 norm
1e2a24a8b Update Github actions to avoid deprecation warnings
285e5f9f4 Do not update brew formulae upon install
73d95b03f Clang-Tidy fixes
51d52c3ea Correct epsilon in CUDA QR test changed by last commit.
546f5337b Updates to CUDA dense linear algebra tests
19ab2c179 BlockRandomAccessMatrix Refactor
a3a062d72 Add a missing header
2b88bedb2 Remove unused variables
9beea728f Fix a bug in CoordinateDescentMinimizer
8e5d83f07 ClangFormat and ClangTidy changes
b15851508 Parallel operations on vectors
2fd81de12 Add build configuration with CUDA on Linux
06bfe6ffa Remove OpenMP and No threading backends.
352b320ab Fixed SuiteSparse 6.0 version parsing
8fd7828e3 ClangTidy fixes
0424615dc Fix PartitionedMatrixView usage in evaluation_benchmark
77a54dd3d Parallel updates to block-diagonal EtE FtF
946fa50de ClangTidy fixes
e4bef9505 Refactor PartitionedMatrixView to cache the partitions
d3201798e Clean up sparse_cholesky_test
addcd342f ClangTidy fixes
c2e7002d2 Remove an unused variable from evaluation_benchmark.cc
fef6d5875 Parallel left products for PartitionedMatrixView
37a3cb384 Update SuiteSparse in MSVC Github workflow
5d53d1ee3 Parallel for with iteration costs and left product
9aa52c6ff Use FindCUDAToolkit for CMake >= 3.17
47e03a6d8 Add const accessor for Problem::Options used by Problem
d6a931009 Clang Tidy Fixes
9364e31ee Fix a regression in SuiteSparse::AnalyzeCholesky
89b3e1f88 Remove unused includes of gflags and gtest
984079003 Fix missing regex dependency for gtest on QNX
6b296f27f Fix missing namespace qualification and docs for Manifold gtest macro
6685e629f AddBlockStructureTranspose to BlockSparseMatrix
699e3f3b3 Fix a link error in evaluation_benchmark.cc
19a3d07f9 Add multiplication benchmarks on BAL data
b221b1294 Format code with clang-format.
ccf32d70c Purge all remaining references to (defunct) LocalParameterization
5f8c406e2 struct ContextImpl -> class ContextImpl
9893c534c Several cleanups.
a78a57472 ClangTidy fixes
b1fe60330 Parallel right products for partitioned view
16668eedf Fix a memory leak in ContextImpl
d129938d5 Third time is the charm
4129b214a More fixes to cuda_dense_cholesky_test.cc
5c01d2573 Remove unused variables from cuda_dense_cholesky_test.cc
5e877ae69 Fix the Bazel build
d89290ffa Fix evalution_benchmark compilability
ae7f456e3 ClangTidy fixes
9438c370f Restore the semantics of TrustRegionMinimizer
c964fce90 Fix bug in cuda_kernels_test
afaad5678 Fix a typo
b7116824b Evaluation benchmark
2b89ce66f Add generalized Euler Angle conversions
8230edc6c ClangTidy fixes
9a2894763 Speed up locking when num_threads = 1.
739f2a25a Parallelize block_jacobi_preconditioner
c0c4f9394 Change implementation of parallel for
fc826c578 CUDA Cleanup
660af905f Fix a bug in TrustRegionMinimizer.
4cd257cf4 Let NumericDiffFirstOrderFunction take a dynamically sized parameter vector
6c27ac6d5 Fix repeated SpMV Benchmark
430a292ac Add a missing include
858b4b89b More ClangTidy fixes
e9e995740 ClangTidy fixes
7f8e930a0 Fix lint errors
f86a3bdbe Unify Block handling across matrix types
5f1946879 clang-formated source
00a05cf70 CUDA SDK Version-based SpMV Selection
de0f74e40 Optimize  the BlockCRSJacobiPreconditioner
ba65ddd31 Improvements to benchmarks
42352e2e2 Added CUDA Jacobi Preconditioner
f802a09ff &foo[0] -> foo.data()
344929647 Add CUDA GPU and Runtime Detection
6085e45be Minor CUDA cleanup.
e15ec89f3 Speed up bundle_adjuster
058a72782 One more ClangTidy fix.
f6f2f0d16 ClangTidy cleanups Also some clang-format cleanups.
9d74c6913 Fix some more errant CATD warnings
a4f744095 Remove an unused variable from compressed_row_sparse_matrix.cc
388c14286 Fix GCC 12.1.1 LTO -Walloc-size-larger-than= warnings
9ec4f7e44 Refactor BlockJacobiPreconditioner
adda97acd Fixed a few more missing CERES_NO_CUDA guards
22aeb3584 Fix a missing string assignment in solver.cc
3b891f767 Insert missing CUDA guards.
829089053 CUDA CGNR, Part 4: CudaCgnrSolver
6ab435d77 Fix a missing CERES_NO_CUDA guard
c560bc2be CUDA CGNR, Part 3: CudaSparseMatrix
c914c7a2b CUDA CGNR, Part 2: CudaVector
3af3dee18 Simplify the implementation to convert from BlockSparseMatrix to CompressedRowSparseMatrix.
242fc0795 Remove unnecessary destructors
737200ac8 Add macos-12 to Github workflow runners
9c968d406 Silence Clang warning
2c78c5f33 Small naming fixups.
2a25d86b0 Integrate schur power series expansion options to bundle adjuster
4642e4b0c Use if constexpr and map SparseMatrix as const
92d837953 Enable usage of schur power series expansion preconditioner.
f1dfac8cd Reduce the number of individual PRNG instances
79e403b15 Expand vcpkg installation instructions
7b0bb0e3f ClangTidy cleanups
c5c2afcc9 Fix solver_test.cc for preconditioners and sparse linear algebra libraries
07d333fb6 Refactor options checking for linear solvers
ba7207b0b A number of small changes.
20e85bbe3 Add power series expansion preconditioner
04899645c LinearOperator::FooMultiply -> LinearOperator::FooMultiplyAndAccumulate
288a3fde6 Add missing virtual destructors to matrix adapters
6483a2b4c Add Sergiu's name to the list of maintainers
1cf49d688 Update FindGlog.cmake to create glog::glog target
1da72ac39 Refactor ConjugateGradientsSolver
f62dccdb3 Fix the Sphere and Line Manifold formulations
3e1cc89f6 A bunch of clang-tidy fixes.
80380538a One more CATD fix
560ef46fb A bunch of minor fixes.
67bae28c1 CUDA CGNR, Part 1: Misc. CLeanup
5d0bca14d Remove ceres/internal/random.h in favor of <random>
d881b5ccf Minor fixes in comments
37516c968 Fix a bug in InnerProductComputer.
d8dad14ee CUDA Cleanup
c9d2ec8a9 Updates to sparse block matrix structures to support new sparse linear solvers.
5fe0bd45a Added MinGW to Windows Github workflow
738c027c1 Fix a logic error in iterative_refiner_test
cb6ad463d Add mixed precision support for CPU based DenseCholesky
df55682ba Fix Eigen error in 2D sphere manifolds Since Eigen does not allow to have a RowMajor column vector (see https://gitlab.com/libeigen/eigen/-/issues/416), the storage order must be set to ColMajor in that case. This fix adds that special case when generating 2D sphere manifolds.
1cf59f61e Set Github workflow NDK path explicitly
68c53bb39 Remove ceres::LocalParameterization
2f660464c Fix build issue with CUDA testing targets when compiling without gflags.
c801192d4 Minor fixes
ce9e902b8 Fix missing CERES_METIS_VERSION
d9a3dfbf2 Add a missing ifdef guard to dense_cholesky_test
5bd43a1fa Speed up DenseSparseMatrix::SquareColumnNorm.
cbc86f651 Fix the build when CUDA is not present
5af8e6449 Update year in solver.h
88e08cfe7 Mixed-precision Iterative Refinement Cholesky With CUDA
290b34ef0 Fix optional SuiteSparse + METIS test-suite names to be unique
d038e2d83 Fix use of NESDIS with SuiteSparse in tests if METIS is not found
027e741a1 Eliminated MinGW warning
4e5ea292b Fixed MSVC 2022 warning
83f6e0853 Fix use of conditional preprocessor checks within a macro in tests
70f1aac31 Fix fmin/fmax() when using Jets with float as their scalar type
5de77f399 Fix reporting of METIS version
11e637667 Fix #ifdef guards around METIS usage in EigenSparse backend and tests
0c88301e6 Provide optional METIS support
f11c25626 Fix fmin/fmax() to use Jet averaging on equality
b90053f1a Revert C++17 usage of std::exclusive_scan
dfce1e128 Link against threading library only if necessary
69eddfb6d Use find module to link against OpenMP
b4803778c Update documentation for linear_solver_ordering_type
2e764df06 Update Cuda memcheck test
443ae9ce2 Update Cuda memcheck test
55b4c3f44 Retain terminal formatting when building docs
786866d9f Generate version string at compile time
5bd83c4ac Unbreak the build with EIGENSPARSE is disabled
2335b5b4b Remove support for CXSparse
fbc2eea16 Nested dissection for ACCELERATE_SPARSE & EIGEN_SPARSE
d87fd551b Fix Ubuntu 20.04 workflow tests
71717f37c Use glog 0.6 release to run Windows Github workflow
66e0adfa7 Fix detection of sphinx-rtd-theme
d09f7e9d5 Enable postordering when computing the sparse factorization.
9b34ecef1 Unbreak the build on MacOS
8ba8fbb17 Remove Solver::Options::use_postordering
30b4d5df3 Fix the ceres.bzl to add missing cc files.
39ec5e8f9 Add Nested Dissection based fill reducing ordering
aa62dd86a Fix a build breakage
41c5fb1e8 Refactor suitesparse.h/cc
12263e283 Make the min. required version of SuiteSparse to be 4.5.6
c8493fc36 Convert internal enums to be class enums.
bb3a40c09 Add Nested Dissection ordering method to SuiteSparse
f1414cb5b Correct spelling in comments and docs.
fd2b0ceed Correct spelling (contiguous, BANS)
464abc198 Run Linux Github workflow on Ubuntu 22.04
caf614a6c Modernize code using c++17 constructs
be618133e Simplify some template metaprograms using fold expressions.
3b0096c1b Add the ability to specify the pivot threshold in Covariance::Options
40c1a7e18 Fix Github workflows
127474360 Ceres Solver now requires C++17
b5f1b7877 clang-format cleanup
32cd1115c Make the code in small_blas_generic.h more compiler friendly.
f68321e7d Update version history
b34280207 Fix MSVC small_blas_test failures
b246991b6 Update the citation instructions in the docs
c0c14abca Fix version history item numbering
d23dbac25 Update Windows install guide
e669c9fc7 Provide citation file
ff57c2e91 Update version history for 2.1.0rc2
ab9436cb9 Workaround MSVC STL deficiency in C++17 mode
97c232857 Update the included gtest to version 1.11.0
4eac7ddd2 Fix Jet lerp test regression
2ffbe126d Fix Jet test failures on ARMv8 with recent Xcode
0d6a0292c Fix unused arguments of Make1stOrderPerturbation
93511bfdc Fix SuiteSparse path and version reporting
bf329b30f Fix link to macOS badge
0133dada2 Add Github workflows
3d3d6ed71 Add missing includes
0a9c0df8a Fix path for cuda-memcheck tests
ee35ef66f ClangFormat cleanup via scripts/all_format.sh
470515985 Add missing includes for config.h
d3612c12c Set CMP0057 policy for IN_LIST operator in FindSuiteSparse.cmake
4bc100c13 Do not define unusable import targets
e91995cce Fix Ubuntu 18.04 shared library build
94af09186 Force C++ linker
a65e73885 Update installation docs
1a377d707 Fix Ubuntu 20.04 SPQR build
817f5a068 Switch to imported SuiteSparse, CXSparse, and METIS targets
b0f32a20d Hide remaining internal symbols
572395098 Add a missing include
b0aef211d Allow to store pointers in ProductManifold
9afe8cc45 Small compile fix to context_impl
284be88ca Allow ProductManifold default construction
f59059fff Bugfix to CUDA workspace handling
779634164 Fix MSVC linker error
7743d2e73 Store ProductManifold instances in a tuple
9f32c42ba Update Travis-CI status badge to .com from .org
eadfead69 Move LineManifold and SphereManifold into their own headers.
f0f8f93bb Fix docs inconsistencies
e40391efa Update version history in preparation for 2.1.0
6a37fbf9b Add static/compile time sizing to EuclideanManifold
4ad787ce1 Fix the bazel build
ae4d95df6 Two small clang-tidy fixes
f0851667b Fix MSVC compilation errors
c8658c899 Modernize more
46b3495a4 Standardize path handling using GNUInstallDirs
98bc3ca17 Fix shared library build due to missing compile features specification
8fe8ebc3a Add final specifier to public classes
84e1696f4 Add final specifier to internal classes.
518970f81 Context should be exported
09ec4997f Cleanup examples
90e58e10f Add missing #include.
15348abe9 Add CUDA based bundle adjustment tests.
57ec9dc92 Do not enforce a specific C++ standard
99698f053 Fix Apple Clang weak symbols warnings
8e0842162 Add support for dense CUDA solvers #3
aff51c907 Revert "Do not enforce a specific C++ standard"
527c3f7da Fixed gflags dependency
d839b7792 Do not enforce a specific C++ standard
f71167c62 Fixed missing include in minilog build
47502b833 Miscellaneous CUDA related changes.
bb2996681 Check CUDA is available in solver.cc
7d2e4152e Add support for dense CUDA solvers #2
f90833f5f Simplify symbol export
c6158e0ab Replace NULL by nullptr
7e4f5a51b Remove blas.h/cc as they are not used anymore.
e0fef6ef0 Add cmake option ENABLE_BITCODE for iOS builds
446487c54 Add <memory> header to all files using std::unique_ptr.
9c5f29d46 Use compiler attributes instead of [[deprecated]]
44039af2c Convert factory functions to return std::unique_ptrs.
708a2a723 Silence LocalParameterization deprecation warnings
de69e657a Fix another missing declaration warning
677711138 Fix segmentation fault in AVX2 builds
0141ca090 Deprecate LocalParameterizations
fdfa5184a Fix missing declaration warning
4742bf386 Add LineManifold.
c14f360e6 Drop trivial special members
ae65219e0 ClangTidy cleanups
a35bd1bf9 Use = default for trivial special members
db67e621e Fix dense_cholesky_test.cc comma handling.
484d3414e Replace virtual keyword by override
2092a720e Fix some nits.
36d6d8690 Add support for dense CUDA solvers #1
af5e48c71 Add SphereManifold.
182cb01c5 Normalize Jet classification and comparison
9dbd28989 Loosen tolerances in dense_qr_test.cc
cab853fd5 Add DenseQR Interface
8ae054ad9 Fix missing declaration warning in autodiff_manifold_test
408af7b1a Move the constructor and destructor for SchurComplementSolver
d51672d1c Move the constructor and destructor for DenseSchurComplementSolver
177b2f99d Use benchmark version 1.6 compatible syntax.
ce9669003 Add const accessor for functor wrapped by auto/numeric-diff objects
9367ec9cc Add a benchmark for dense linear solvers.
7d6524daf Support fma Jet
a0d81ad63 Fix a bug in AutoDiffManifold
5a99e42e1 ClangTidy fixes
40fb41355 Update .gitignore
e6e6ae087 Unbreak the bazel build
0572efc57 Fix a compilation warning in autodiff_manifold.h
475db73d0 Fix build breakage when LAPACK support is disabled.
095c9197f Fix iterative_refiner_test.cc
6d06e9b98 Add DenseCholesky
77c0c4d09 Migrate examples to use Manifolds
19eef54fc Rename Quaternion to QuaternionManifold. Also rename EigenQuaternion to EigenQuaternionManiold.
ca6d841c2 Add AutoDiffManifold
97d7e0737 Move the manifold testing matchers to manifold_test_utils.h
16436b34b Fix some more clang-tidy suggestions.
dcdefc216 Fix a bunch of clang-tidy suggestions.
d8a1b69ab Remove an unused variable from gradient_checker.cc
611b46b54 Remove the use of CHECK_NOTNULL.
125a0e9be LocalParameterization -> Manifold #1
00bfbae11 Add missing algorithm header to manifold.cc
1d5aff059 Refactor jet_tests.cc
fbd693091 Fix two unused variable warnings.
c0cb42e5f Add Problem::HasParameterization
7e2f9d9d4 Add EigenQuaternion manifold and tests for it.
b81a8bbb7 Add the Quaternion manifold and tests.
4a01dcb88 Add more invariants and documentation to manifold_test.cc
bdd80fcce Improve Manifold testing
ce1537030 Fixed missing headers in manifold.h
23b204d7e LocalParameterization -> Manifold #1
c2fab6502 Fix docs of supported sparse backends for mixed_precision_solves option
8cb441c49 Fix missing declaration warnings in GCC
d2b7f337c Remove split.cc from the bazel sources.
206061a6b Use standard c++ types in jet_test.cc
1f374a9fe Support promotion in comparison between Jet and scalars
06e68dbc5 Avoid midpoint overflow in the differential
276d24c73 Fix C++20 compilation
b1391e062 Support midpoint Jet
8426526df Support lerp Jet
57c279689 support 3-argument hypot jet
123fba61c Eigen::MappedSparseMatrix -> Eigen::Map<Eigen::SparseMatrix>
3f950c66d reworked copysign tests
48cb54d1b fix fmin and fmax NaN handling
552a4e517 support log10 jet
4e49c5422 reuse expm1 result for differential
8d3e64dd5 Use modern-style Eigen3 CMake variables
2fba61434 support log1p and expm1 jet
a668cabbc support norm jet
a3a4b6d77 support copysign jet
b75dac169 fix abs jet test comment
034bf566f use copysign for abs jet
31008453f Add example for BiCubicInterpolator
7ef4a1221 Add a section on implicit and inverse function theorems
686428f5c Move the further reading section to bibliography
e47d87fdd Add a note about Trigg's correction
d2852518d Fix the docs for Problem::RemoveResidualBlock & Problem::RemoveParameterBlock
06e02a173 Delete unused files split.h/cc
ac7268da5 Fix an 80cols issue in covariance_impl.cc
17dccef91 Add NumericDiffFirstOrderFunction
03d64141a Fix a incorrect check in reorder_program.cc
884111913 Two changes to TinySolver
4dff3ea2c Fix a number of typos in rotation.h
98719ced4 Fix a type in interfacing_with_autodiff.html
0299ce944 Update conf.py to be compatible with Sphinx 4.1.2
2a2b9bd6f Fix a bug in covariance_impl.cc
27fade7b8 Fix a bug in system_test.cc
d4eb83ee5 Fix the Jacobian in trust_region_minimizer_test.cc
5f6071a1c Fix a bug in local_parameterization_test.cc
b2e732b1e Fix errors in comments from William Gandler.
42f1d6717 Add accessors to GradientProblem
aefd37b18 Refactor small_blas_gemm_benchmark
dc20db303 [docs] Fix `IterationSummary` section. Add missing `IterationCallback`
90ba7d1ef [docs] Fix typos
c3129c3d4 Fix tests not executing
7de561e8e Fix dependency check for building documentation
4fbe218f2 Refactor small_blas_test
3fdde8ede Remove an errant double link.
20ad431f7 Fixing a typo in the version history
0c85c4092 Revert "Reduce copies involved in Jet operations"
3a02d5aaf Fix typo in LossFunctionWrapper sample code
7b2c223be Add fmax/fmin overloads for scalars
c036c7819 Reduce copies involved in Jet operations
51945e061 Introduce benchmark for Jet operations
ec4f2995b Do not check MaxNumThreadsAvailable if the thread number is set to 1.
98f639f54 Add a macro CERES_GET_FLAG.
766f2cab5 Reduce log spam in covariance_impl.cc.
941ea1347 Fix FindTBB version detection with TBB >= 2021.1.1
323c350a6 fix Eigen3_VERSION
2b32b3212 Revert "Group specializations into groups of four"
313caf1ae Allow Unity build.
4ba244cdb Group specializations into groups of four
d77a8100a Make miniglog's InitGoogleLogging argument const.
863724994 Use portable expression for constant 2/sqrt(pi)
97873ea65 Add some missing includes for glog/logging.h
d15b1bcd3 Increase tolerance in small_blas_test.cc
17cf01831 Hide 'format not a string literal' error in examples
64029909b Fix -Wno-maybe-uninitialized error
21294123d Fix nonnull arg compared to NULL error.
1dd417410 Fix -Wno-format-nonliteral
6c106bf51 Fix -Wmissing-field-initializers error
c48a32792 Use cc_binary includes so examples build as external repo
e0e14a5cd Fix errors found by -Werror
e84cf10e1 Fix an errant double in TinySolver.
66b4c33e8 updated unit quaternion rotation
d45ec47b5 Fix a typo in schur_eliminator.h

Change-Id: I60db062e44d051d50dbb3a145eec2f74d5190481
git-subtree-dir: third_party/ceres
git-subtree-split: 7f9cc571b03632f1df93ea35725a1f5dfffe2c72
Signed-off-by: Austin Schuh <austin.linux@gmail.com>
diff --git a/examples/BUILD b/examples/BUILD
index 90daf21..4f68eda 100644
--- a/examples/BUILD
+++ b/examples/BUILD
@@ -31,6 +31,8 @@
 EXAMPLE_COPTS = [
     # Needed to silence GFlags complaints.
     "-Wno-sign-compare",
+    # Needed to put fscanf in a function.
+    "-Wno-format-nonliteral",
 ]
 
 EXAMPLE_DEPS = [
@@ -45,7 +47,6 @@
         "bal_problem.cc",
         "bal_problem.h",
         "bundle_adjuster.cc",
-        "random.h",
         "snavely_reprojection_error.h",
     ],
     copts = EXAMPLE_COPTS,
@@ -67,26 +68,24 @@
 cc_binary(
     name = "robot_pose_mle",
     srcs = [
-        "random.h",
         "robot_pose_mle.cc",
     ],
     copts = EXAMPLE_COPTS,
     deps = EXAMPLE_DEPS,
 )
 
-SLAM_COPTS = EXAMPLE_COPTS + ["-Iexamples/slam"]
-
 cc_binary(
     name = "pose_graph_2d",
     srcs = [
         "slam/common/read_g2o.h",
-        "slam/pose_graph_2d/angle_local_parameterization.h",
+        "slam/pose_graph_2d/angle_manifold.h",
         "slam/pose_graph_2d/normalize_angle.h",
         "slam/pose_graph_2d/pose_graph_2d.cc",
         "slam/pose_graph_2d/pose_graph_2d_error_term.h",
         "slam/pose_graph_2d/types.h",
     ],
-    copts = SLAM_COPTS,
+    copts = EXAMPLE_COPTS,
+    includes = ["slam"],
     deps = EXAMPLE_DEPS,
 )
 
@@ -98,7 +97,8 @@
         "slam/pose_graph_3d/pose_graph_3d_error_term.h",
         "slam/pose_graph_3d/types.h",
     ],
-    copts = SLAM_COPTS,
+    copts = EXAMPLE_COPTS,
+    includes = ["slam"],
     deps = EXAMPLE_DEPS,
 )
 
diff --git a/examples/CMakeLists.txt b/examples/CMakeLists.txt
index 7f9b117..8af2077 100644
--- a/examples/CMakeLists.txt
+++ b/examples/CMakeLists.txt
@@ -1,5 +1,5 @@
 # Ceres Solver - A fast non-linear least squares minimizer
-# Copyright 2015 Google Inc. All rights reserved.
+# Copyright 2022 Google Inc. All rights reserved.
 # http://ceres-solver.org/
 #
 # Redistribution and use in source and binary forms, with or without
@@ -28,86 +28,96 @@
 #
 # Author: keir@google.com (Keir Mierle)
 
-# Only Ceres itself should be compiled with CERES_BUILDING_SHARED_LIBRARY
-# defined, any users of Ceres will have CERES_USING_SHARED_LIBRARY defined
-# for them in Ceres' config.h if appropriate.
-if (BUILD_SHARED_LIBS)
-  remove_definitions(-DCERES_BUILDING_SHARED_LIBRARY)
-endif()
-
 add_executable(helloworld helloworld.cc)
-target_link_libraries(helloworld Ceres::ceres)
+target_link_libraries(helloworld PRIVATE Ceres::ceres)
 
 add_executable(helloworld_numeric_diff helloworld_numeric_diff.cc)
-target_link_libraries(helloworld_numeric_diff Ceres::ceres)
+target_link_libraries(helloworld_numeric_diff PRIVATE Ceres::ceres)
 
 add_executable(helloworld_analytic_diff helloworld_analytic_diff.cc)
-target_link_libraries(helloworld_analytic_diff Ceres::ceres)
+target_link_libraries(helloworld_analytic_diff PRIVATE Ceres::ceres)
 
 add_executable(curve_fitting curve_fitting.cc)
-target_link_libraries(curve_fitting Ceres::ceres)
+target_link_libraries(curve_fitting PRIVATE Ceres::ceres)
 
 add_executable(rosenbrock rosenbrock.cc)
-target_link_libraries(rosenbrock Ceres::ceres)
+target_link_libraries(rosenbrock PRIVATE Ceres::ceres)
+
+add_executable(rosenbrock_analytic_diff rosenbrock_analytic_diff.cc)
+target_link_libraries(rosenbrock_analytic_diff PRIVATE Ceres::ceres)
+
+add_executable(rosenbrock_numeric_diff rosenbrock_numeric_diff.cc)
+target_link_libraries(rosenbrock_numeric_diff PRIVATE Ceres::ceres)
 
 add_executable(curve_fitting_c curve_fitting.c)
-target_link_libraries(curve_fitting_c Ceres::ceres)
-# Force CMake to link curve_fitting_c using the C linker.
-set_target_properties(curve_fitting_c PROPERTIES LINKER_LANGUAGE C)
+target_link_libraries(curve_fitting_c PRIVATE Ceres::ceres)
+# Force CMake to link curve_fitting_c using the C++ linker.
+set_target_properties(curve_fitting_c PROPERTIES LINKER_LANGUAGE CXX)
 # As this is a C file #including <math.h> we have to explicitly add the math
 # library (libm). Although some compilers (dependent upon options) will accept
 # the indirect link to libm via Ceres, at least GCC 4.8 on pure Debian won't.
-if (NOT MSVC)
-  target_link_libraries(curve_fitting_c m)
-endif (NOT MSVC)
+if (HAVE_LIBM)
+  target_link_libraries(curve_fitting_c PRIVATE m)
+endif (HAVE_LIBM)
 
 add_executable(ellipse_approximation ellipse_approximation.cc)
-target_link_libraries(ellipse_approximation Ceres::ceres)
+target_link_libraries(ellipse_approximation PRIVATE Ceres::ceres)
 
 add_executable(robust_curve_fitting robust_curve_fitting.cc)
-target_link_libraries(robust_curve_fitting Ceres::ceres)
+target_link_libraries(robust_curve_fitting PRIVATE Ceres::ceres)
 
 add_executable(simple_bundle_adjuster simple_bundle_adjuster.cc)
-target_link_libraries(simple_bundle_adjuster Ceres::ceres)
+target_link_libraries(simple_bundle_adjuster PRIVATE Ceres::ceres)
+
+add_executable(bicubic_interpolation bicubic_interpolation.cc)
+target_link_libraries(bicubic_interpolation PRIVATE Ceres::ceres)
+
+add_executable(bicubic_interpolation_analytic bicubic_interpolation_analytic.cc)
+target_link_libraries(bicubic_interpolation_analytic PRIVATE Ceres::ceres)
+
+add_executable(iteration_callback_example iteration_callback_example.cc)
+target_link_libraries(iteration_callback_example PRIVATE Ceres::ceres)
+
+add_executable(evaluation_callback_example evaluation_callback_example.cc)
+target_link_libraries(evaluation_callback_example PRIVATE Ceres::ceres)
 
 if (GFLAGS)
   add_executable(powell powell.cc)
-  target_link_libraries(powell Ceres::ceres gflags)
+  target_link_libraries(powell PRIVATE Ceres::ceres gflags)
 
   add_executable(nist nist.cc)
-  target_link_libraries(nist Ceres::ceres gflags)
-  if (MSVC)
-    target_compile_options(nist PRIVATE "/bigobj")
+  target_link_libraries(nist PRIVATE Ceres::ceres gflags)
+  if (HAVE_BIGOBJ)
+    target_compile_options(nist PRIVATE /bigobj)
   endif()
 
   add_executable(more_garbow_hillstrom more_garbow_hillstrom.cc)
-  target_link_libraries(more_garbow_hillstrom Ceres::ceres gflags)
+  target_link_libraries(more_garbow_hillstrom PRIVATE Ceres::ceres gflags)
 
   add_executable(circle_fit circle_fit.cc)
-  target_link_libraries(circle_fit Ceres::ceres gflags)
+  target_link_libraries(circle_fit PRIVATE Ceres::ceres gflags)
 
   add_executable(bundle_adjuster
                  bundle_adjuster.cc
                  bal_problem.cc)
-  target_link_libraries(bundle_adjuster Ceres::ceres gflags)
+  target_link_libraries(bundle_adjuster PRIVATE Ceres::ceres gflags)
 
   add_executable(libmv_bundle_adjuster
                  libmv_bundle_adjuster.cc)
-  target_link_libraries(libmv_bundle_adjuster Ceres::ceres gflags)
+  target_link_libraries(libmv_bundle_adjuster PRIVATE Ceres::ceres gflags)
 
   add_executable(libmv_homography
                  libmv_homography.cc)
-  target_link_libraries(libmv_homography Ceres::ceres gflags)
+  target_link_libraries(libmv_homography PRIVATE Ceres::ceres gflags)
 
   add_executable(denoising
                  denoising.cc
                  fields_of_experts.cc)
-  target_link_libraries(denoising Ceres::ceres gflags)
+  target_link_libraries(denoising PRIVATE Ceres::ceres gflags)
 
   add_executable(robot_pose_mle
                  robot_pose_mle.cc)
-  target_link_libraries(robot_pose_mle Ceres::ceres gflags)
-
+  target_link_libraries(robot_pose_mle PRIVATE Ceres::ceres gflags)
 endif (GFLAGS)
 
 add_subdirectory(sampled_function)
diff --git a/examples/Makefile.example b/examples/Makefile.example
deleted file mode 100644
index f2b0dc0..0000000
--- a/examples/Makefile.example
+++ /dev/null
@@ -1,82 +0,0 @@
-# Ceres Solver - A fast non-linear least squares minimizer
-# Copyright 2015 Google Inc. All rights reserved.
-# http://ceres-solver.org/
-#
-# Redistribution and use in source and binary forms, with or without
-# modification, are permitted provided that the following conditions are met:
-#
-# * Redistributions of source code must retain the above copyright notice,
-#   this list of conditions and the following disclaimer.
-# * Redistributions in binary form must reproduce the above copyright notice,
-#   this list of conditions and the following disclaimer in the documentation
-#   and/or other materials provided with the distribution.
-# * Neither the name of Google Inc. nor the names of its contributors may be
-#   used to endorse or promote products derived from this software without
-#   specific prior written permission.
-#
-# THIS SOFTWARE IS PROVIDED BY THE COPYRIGHT HOLDERS AND CONTRIBUTORS "AS IS"
-# AND ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED TO, THE
-# IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR PURPOSE
-# ARE DISCLAIMED. IN NO EVENT SHALL THE COPYRIGHT OWNER OR CONTRIBUTORS BE
-# LIABLE FOR ANY DIRECT, INDIRECT, INCIDENTAL, SPECIAL, EXEMPLARY, OR
-# CONSEQUENTIAL DAMAGES (INCLUDING, BUT NOT LIMITED TO, PROCUREMENT OF
-# SUBSTITUTE GOODS OR SERVICES; LOSS OF USE, DATA, OR PROFITS; OR BUSINESS
-# INTERRUPTION) HOWEVER CAUSED AND ON ANY THEORY OF LIABILITY, WHETHER IN
-# CONTRACT, STRICT LIABILITY, OR TORT (INCLUDING NEGLIGENCE OR OTHERWISE)
-# ARISING IN ANY WAY OUT OF THE USE OF THIS SOFTWARE, EVEN IF ADVISED OF THE
-# POSSIBILITY OF SUCH DAMAGE.
-#
-# Author: keir@google.com (Keir Mierle)
-#
-# This is an example Makefile for using Ceres. In practice, the Ceres authors
-# suggest using CMake instead, but if Make is needed for some reason, this
-# example serves to make it easy to do so.
-
-# This should point to place where you unpacked or cloned Ceres.
-CERES_SRC_DIR := /home/keir/wrk/ceres-extra
-
-# This should point to the place where you built Ceres. If you got Ceres by
-# installing it, then this will likely be /usr/local/lib.
-CERES_BIN_DIR := /home/keir/wrk/ceres-extra-bin
-
-# The place you unpacked or cloned Eigen. If Eigen was installed from packages,
-# this will likely be /usr/local/include.
-EIGEN_SRC_DIR := /home/keir/src/eigen-3.0.5
-
-INCLUDES := -I$(CERES_SRC_DIR)/include \
-            -I$(EIGEN_SRC_DIR)
-
-CERES_LIBRARY := -lceres
-CERES_LIBRARY_PATH := -L$(CERES_BIN_DIR)/lib
-CERES_LIBRARY_DEPENDENCIES = -lgflags -lglog
-
-# If Ceres was built with Suitesparse:
-CERES_LIBRARY_DEPENDENCIES += -llapack -lcamd -lamd -lccolamd -lcolamd -lcholmod
-
-# If Ceres was built with CXSparse:
-CERES_LIBRARY_DEPENDENCIES += -lcxsparse
-
-# If Ceres was built with OpenMP:
-CERES_LIBRARY_DEPENDENCIES += -fopenmp -lpthread -lgomp -lm
-
-# The set of object files for your application.
-APPLICATION_OBJS := simple_bundle_adjuster.o
-
-all: simple_bundle_adjuster
-
-simple_bundle_adjuster: $(APPLICATION_OBJS)
-	g++ \
-		$(APPLICATION_OBJS) \
-		$(CERES_LIBRARY_PATH) \
-		$(CERES_LIBRARY) \
-		$(CERES_LIBRARY_DEPENDENCIES) \
-		-o simple_bundle_adjuster
-
-# Disabling debug asserts via -DNDEBUG helps make Eigen faster, at the cost of
-# not getting handy assert failures when there are issues in your code.
-CFLAGS := -O2 -DNDEBUG
-
-# If you have files ending in .cpp instead of .cc, fix the next line
-# appropriately.
-%.o: %.cc $(DEPS)
-	g++ -c -o $@ $< $(CFLAGS) $(INCLUDES)
diff --git a/examples/bal_problem.cc b/examples/bal_problem.cc
index ceac89a..ccf7449 100644
--- a/examples/bal_problem.cc
+++ b/examples/bal_problem.cc
@@ -1,5 +1,5 @@
 // Ceres Solver - A fast non-linear least squares minimizer
-// Copyright 2015 Google Inc. All rights reserved.
+// Copyright 2023 Google Inc. All rights reserved.
 // http://ceres-solver.org/
 //
 // Redistribution and use in source and binary forms, with or without
@@ -30,22 +30,22 @@
 
 #include "bal_problem.h"
 
+#include <algorithm>
 #include <cstdio>
-#include <cstdlib>
 #include <fstream>
+#include <functional>
+#include <random>
 #include <string>
 #include <vector>
 
 #include "Eigen/Core"
 #include "ceres/rotation.h"
 #include "glog/logging.h"
-#include "random.h"
 
-namespace ceres {
-namespace examples {
+namespace ceres::examples {
 namespace {
-typedef Eigen::Map<Eigen::VectorXd> VectorRef;
-typedef Eigen::Map<const Eigen::VectorXd> ConstVectorRef;
+using VectorRef = Eigen::Map<Eigen::VectorXd>;
+using ConstVectorRef = Eigen::Map<const Eigen::VectorXd>;
 
 template <typename T>
 void FscanfOrDie(FILE* fptr, const char* format, T* value) {
@@ -55,15 +55,14 @@
   }
 }
 
-void PerturbPoint3(const double sigma, double* point) {
+void PerturbPoint3(std::function<double()> dist, double* point) {
   for (int i = 0; i < 3; ++i) {
-    point[i] += RandNormal() * sigma;
+    point[i] += dist();
   }
 }
 
 double Median(std::vector<double>* data) {
-  int n = data->size();
-  std::vector<double>::iterator mid_point = data->begin() + n / 2;
+  auto mid_point = data->begin() + data->size() / 2;
   std::nth_element(data->begin(), mid_point, data->end());
   return *mid_point;
 }
@@ -73,12 +72,12 @@
 BALProblem::BALProblem(const std::string& filename, bool use_quaternions) {
   FILE* fptr = fopen(filename.c_str(), "r");
 
-  if (fptr == NULL) {
+  if (fptr == nullptr) {
     LOG(FATAL) << "Error: unable to open file " << filename;
     return;
   };
 
-  // This wil die horribly on invalid files. Them's the breaks.
+  // This will die horribly on invalid files. Them's the breaks.
   FscanfOrDie(fptr, "%d", &num_cameras_);
   FscanfOrDie(fptr, "%d", &num_points_);
   FscanfOrDie(fptr, "%d", &num_observations_);
@@ -111,7 +110,7 @@
   if (use_quaternions) {
     // Switch the angle-axis rotations to quaternions.
     num_parameters_ = 10 * num_cameras_ + 3 * num_points_;
-    double* quaternion_parameters = new double[num_parameters_];
+    auto* quaternion_parameters = new double[num_parameters_];
     double* original_cursor = parameters_;
     double* quaternion_cursor = quaternion_parameters;
     for (int i = 0; i < num_cameras_; ++i) {
@@ -137,7 +136,7 @@
 void BALProblem::WriteToFile(const std::string& filename) const {
   FILE* fptr = fopen(filename.c_str(), "w");
 
-  if (fptr == NULL) {
+  if (fptr == nullptr) {
     LOG(FATAL) << "Error: unable to open file " << filename;
     return;
   };
@@ -161,8 +160,8 @@
     } else {
       memcpy(angleaxis, parameters_ + 9 * i, 9 * sizeof(double));
     }
-    for (int j = 0; j < 9; ++j) {
-      fprintf(fptr, "%.16g\n", angleaxis[j]);
+    for (double coeff : angleaxis) {
+      fprintf(fptr, "%.16g\n", coeff);
     }
   }
 
@@ -297,14 +296,20 @@
   CHECK_GE(point_sigma, 0.0);
   CHECK_GE(rotation_sigma, 0.0);
   CHECK_GE(translation_sigma, 0.0);
-
+  std::mt19937 prng;
+  std::normal_distribution<double> point_noise_distribution(0.0, point_sigma);
   double* points = mutable_points();
   if (point_sigma > 0) {
     for (int i = 0; i < num_points_; ++i) {
-      PerturbPoint3(point_sigma, points + 3 * i);
+      PerturbPoint3(std::bind(point_noise_distribution, std::ref(prng)),
+                    points + 3 * i);
     }
   }
 
+  std::normal_distribution<double> rotation_noise_distribution(0.0,
+                                                               point_sigma);
+  std::normal_distribution<double> translation_noise_distribution(
+      0.0, translation_sigma);
   for (int i = 0; i < num_cameras_; ++i) {
     double* camera = mutable_cameras() + camera_block_size() * i;
 
@@ -314,12 +319,14 @@
     // representation.
     CameraToAngleAxisAndCenter(camera, angle_axis, center);
     if (rotation_sigma > 0.0) {
-      PerturbPoint3(rotation_sigma, angle_axis);
+      PerturbPoint3(std::bind(rotation_noise_distribution, std::ref(prng)),
+                    angle_axis);
     }
     AngleAxisAndCenterToCamera(angle_axis, center, camera);
 
     if (translation_sigma > 0.0) {
-      PerturbPoint3(translation_sigma, camera + camera_block_size() - 6);
+      PerturbPoint3(std::bind(translation_noise_distribution, std::ref(prng)),
+                    camera + camera_block_size() - 6);
     }
   }
 }
@@ -331,5 +338,4 @@
   delete[] parameters_;
 }
 
-}  // namespace examples
-}  // namespace ceres
+}  // namespace ceres::examples
diff --git a/examples/bal_problem.h b/examples/bal_problem.h
index e6d4ace..93effc3 100644
--- a/examples/bal_problem.h
+++ b/examples/bal_problem.h
@@ -1,5 +1,5 @@
 // Ceres Solver - A fast non-linear least squares minimizer
-// Copyright 2015 Google Inc. All rights reserved.
+// Copyright 2023 Google Inc. All rights reserved.
 // http://ceres-solver.org/
 //
 // Redistribution and use in source and binary forms, with or without
@@ -39,8 +39,7 @@
 
 #include <string>
 
-namespace ceres {
-namespace examples {
+namespace ceres::examples {
 
 class BALProblem {
  public:
@@ -105,7 +104,6 @@
   double* parameters_;
 };
 
-}  // namespace examples
-}  // namespace ceres
+}  // namespace ceres::examples
 
 #endif  // CERES_EXAMPLES_BAL_PROBLEM_H_
diff --git a/examples/bicubic_interpolation.cc b/examples/bicubic_interpolation.cc
new file mode 100644
index 0000000..21b3c7e
--- /dev/null
+++ b/examples/bicubic_interpolation.cc
@@ -0,0 +1,155 @@
+// Ceres Solver - A fast non-linear least squares minimizer
+// Copyright 2023 Google Inc. All rights reserved.
+// http://ceres-solver.org/
+//
+// Redistribution and use in source and binary forms, with or without
+// modification, are permitted provided that the following conditions are met:
+//
+// * Redistributions of source code must retain the above copyright notice,
+//   this list of conditions and the following disclaimer.
+// * Redistributions in binary form must reproduce the above copyright notice,
+//   this list of conditions and the following disclaimer in the documentation
+//   and/or other materials provided with the distribution.
+// * Neither the name of Google Inc. nor the names of its contributors may be
+//   used to endorse or promote products derived from this software without
+//   specific prior written permission.
+//
+// THIS SOFTWARE IS PROVIDED BY THE COPYRIGHT HOLDERS AND CONTRIBUTORS "AS IS"
+// AND ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED TO, THE
+// IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR PURPOSE
+// ARE DISCLAIMED. IN NO EVENT SHALL THE COPYRIGHT OWNER OR CONTRIBUTORS BE
+// LIABLE FOR ANY DIRECT, INDIRECT, INCIDENTAL, SPECIAL, EXEMPLARY, OR
+// CONSEQUENTIAL DAMAGES (INCLUDING, BUT NOT LIMITED TO, PROCUREMENT OF
+// SUBSTITUTE GOODS OR SERVICES; LOSS OF USE, DATA, OR PROFITS; OR BUSINESS
+// INTERRUPTION) HOWEVER CAUSED AND ON ANY THEORY OF LIABILITY, WHETHER IN
+// CONTRACT, STRICT LIABILITY, OR TORT (INCLUDING NEGLIGENCE OR OTHERWISE)
+// ARISING IN ANY WAY OUT OF THE USE OF THIS SOFTWARE, EVEN IF ADVISED OF THE
+// POSSIBILITY OF SUCH DAMAGE.
+//
+// Bicubic interpolation with automatic differentiation
+//
+// We will use estimation of 2d shift as a sample problem for bicubic
+// interpolation.
+//
+// Let us define f(x, y) = x * x - y * x + y * y
+// And optimize cost function sum_i [f(x_i + s_x, y_i + s_y) - v_i]^2
+//
+// Bicubic interpolation of f(x, y) will be exact, thus we can expect close to
+// perfect convergence
+
+#include <utility>
+
+#include "ceres/ceres.h"
+#include "ceres/cubic_interpolation.h"
+#include "glog/logging.h"
+
+using Grid = ceres::Grid2D<double>;
+using Interpolator = ceres::BiCubicInterpolator<Grid>;
+
+// Cost-function using autodiff interface of BiCubicInterpolator
+struct AutoDiffBiCubicCost {
+  EIGEN_MAKE_ALIGNED_OPERATOR_NEW;
+
+  template <typename T>
+  bool operator()(const T* s, T* residual) const {
+    using Vector2T = Eigen::Matrix<T, 2, 1>;
+    Eigen::Map<const Vector2T> shift(s);
+
+    const Vector2T point = point_ + shift;
+
+    T v;
+    interpolator_.Evaluate(point.y(), point.x(), &v);
+
+    *residual = v - value_;
+    return true;
+  }
+
+  AutoDiffBiCubicCost(const Interpolator& interpolator,
+                      Eigen::Vector2d point,
+                      double value)
+      : point_(std::move(point)), value_(value), interpolator_(interpolator) {}
+
+  static ceres::CostFunction* Create(const Interpolator& interpolator,
+                                     const Eigen::Vector2d& point,
+                                     double value) {
+    return new ceres::AutoDiffCostFunction<AutoDiffBiCubicCost, 1, 2>(
+        interpolator, point, value);
+  }
+
+  const Eigen::Vector2d point_;
+  const double value_;
+  const Interpolator& interpolator_;
+};
+
+// Function for input data generation
+static double f(const double& x, const double& y) {
+  return x * x - y * x + y * y;
+}
+
+int main(int argc, char** argv) {
+  google::InitGoogleLogging(argv[0]);
+  // Problem sizes
+  const int kGridRowsHalf = 9;
+  const int kGridColsHalf = 11;
+  const int kGridRows = 2 * kGridRowsHalf + 1;
+  const int kGridCols = 2 * kGridColsHalf + 1;
+  const int kPoints = 4;
+
+  const Eigen::Vector2d shift(1.234, 2.345);
+  const std::array<Eigen::Vector2d, kPoints> points = {
+      Eigen::Vector2d{-2., -3.},
+      Eigen::Vector2d{-2., 3.},
+      Eigen::Vector2d{2., 3.},
+      Eigen::Vector2d{2., -3.}};
+
+  // Data is a row-major array of kGridRows x kGridCols values of function
+  // f(x, y) on the grid, with x in {-kGridColsHalf, ..., +kGridColsHalf},
+  // and y in {-kGridRowsHalf, ..., +kGridRowsHalf}
+  double data[kGridRows * kGridCols];
+  for (int i = 0; i < kGridRows; ++i) {
+    for (int j = 0; j < kGridCols; ++j) {
+      // Using row-major order
+      int index = i * kGridCols + j;
+      double y = i - kGridRowsHalf;
+      double x = j - kGridColsHalf;
+
+      data[index] = f(x, y);
+    }
+  }
+  const Grid grid(data,
+                  -kGridRowsHalf,
+                  kGridRowsHalf + 1,
+                  -kGridColsHalf,
+                  kGridColsHalf + 1);
+  const Interpolator interpolator(grid);
+
+  Eigen::Vector2d shift_estimate(3.1415, 1.337);
+
+  ceres::Problem problem;
+  problem.AddParameterBlock(shift_estimate.data(), 2);
+
+  for (const auto& p : points) {
+    const Eigen::Vector2d shifted = p + shift;
+
+    const double v = f(shifted.x(), shifted.y());
+    problem.AddResidualBlock(AutoDiffBiCubicCost::Create(interpolator, p, v),
+                             nullptr,
+                             shift_estimate.data());
+  }
+
+  ceres::Solver::Options options;
+  options.minimizer_progress_to_stdout = true;
+
+  ceres::Solver::Summary summary;
+  ceres::Solve(options, &problem, &summary);
+  std::cout << summary.BriefReport() << '\n';
+
+  std::cout << "Bicubic interpolation with automatic derivatives:\n";
+  std::cout << "Estimated shift: " << shift_estimate.transpose()
+            << ", ground-truth: " << shift.transpose()
+            << " (error: " << (shift_estimate - shift).transpose() << ")"
+            << std::endl;
+
+  CHECK_LT((shift_estimate - shift).norm(), 1e-9);
+  return 0;
+}
diff --git a/examples/bicubic_interpolation_analytic.cc b/examples/bicubic_interpolation_analytic.cc
new file mode 100644
index 0000000..4b79d56
--- /dev/null
+++ b/examples/bicubic_interpolation_analytic.cc
@@ -0,0 +1,166 @@
+// Ceres Solver - A fast non-linear least squares minimizer
+// Copyright 2023 Google Inc. All rights reserved.
+// http://ceres-solver.org/
+//
+// Redistribution and use in source and binary forms, with or without
+// modification, are permitted provided that the following conditions are met:
+//
+// * Redistributions of source code must retain the above copyright notice,
+//   this list of conditions and the following disclaimer.
+// * Redistributions in binary form must reproduce the above copyright notice,
+//   this list of conditions and the following disclaimer in the documentation
+//   and/or other materials provided with the distribution.
+// * Neither the name of Google Inc. nor the names of its contributors may be
+//   used to endorse or promote products derived from this software without
+//   specific prior written permission.
+//
+// THIS SOFTWARE IS PROVIDED BY THE COPYRIGHT HOLDERS AND CONTRIBUTORS "AS IS"
+// AND ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED TO, THE
+// IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR PURPOSE
+// ARE DISCLAIMED. IN NO EVENT SHALL THE COPYRIGHT OWNER OR CONTRIBUTORS BE
+// LIABLE FOR ANY DIRECT, INDIRECT, INCIDENTAL, SPECIAL, EXEMPLARY, OR
+// CONSEQUENTIAL DAMAGES (INCLUDING, BUT NOT LIMITED TO, PROCUREMENT OF
+// SUBSTITUTE GOODS OR SERVICES; LOSS OF USE, DATA, OR PROFITS; OR BUSINESS
+// INTERRUPTION) HOWEVER CAUSED AND ON ANY THEORY OF LIABILITY, WHETHER IN
+// CONTRACT, STRICT LIABILITY, OR TORT (INCLUDING NEGLIGENCE OR OTHERWISE)
+// ARISING IN ANY WAY OUT OF THE USE OF THIS SOFTWARE, EVEN IF ADVISED OF THE
+// POSSIBILITY OF SUCH DAMAGE.
+//
+// Bicubic interpolation with analytic differentiation
+//
+// We will use estimation of 2d shift as a sample problem for bicubic
+// interpolation.
+//
+// Let us define f(x, y) = x * x - y * x + y * y
+// And optimize cost function sum_i [f(x_i + s_x, y_i + s_y) - v_i]^2
+//
+// Bicubic interpolation of f(x, y) will be exact, thus we can expect close to
+// perfect convergence
+
+#include <utility>
+
+#include "ceres/ceres.h"
+#include "ceres/cubic_interpolation.h"
+#include "glog/logging.h"
+
+using Grid = ceres::Grid2D<double>;
+using Interpolator = ceres::BiCubicInterpolator<Grid>;
+
+// Cost-function using analytic interface of BiCubicInterpolator
+struct AnalyticBiCubicCost : public ceres::CostFunction {
+  EIGEN_MAKE_ALIGNED_OPERATOR_NEW;
+
+  bool Evaluate(double const* const* parameters,
+                double* residuals,
+                double** jacobians) const override {
+    Eigen::Map<const Eigen::Vector2d> shift(parameters[0]);
+
+    const Eigen::Vector2d point = point_ + shift;
+
+    double* f = residuals;
+    double* dfdr = nullptr;
+    double* dfdc = nullptr;
+    if (jacobians && jacobians[0]) {
+      dfdc = jacobians[0];
+      dfdr = dfdc + 1;
+    }
+
+    interpolator_.Evaluate(point.y(), point.x(), f, dfdr, dfdc);
+
+    if (residuals) {
+      *f -= value_;
+    }
+    return true;
+  }
+
+  AnalyticBiCubicCost(const Interpolator& interpolator,
+                      Eigen::Vector2d point,
+                      double value)
+      : point_(std::move(point)), value_(value), interpolator_(interpolator) {
+    set_num_residuals(1);
+    *mutable_parameter_block_sizes() = {2};
+  }
+
+  static ceres::CostFunction* Create(const Interpolator& interpolator,
+                                     const Eigen::Vector2d& point,
+                                     double value) {
+    return new AnalyticBiCubicCost(interpolator, point, value);
+  }
+
+  const Eigen::Vector2d point_;
+  const double value_;
+  const Interpolator& interpolator_;
+};
+
+// Function for input data generation
+static double f(const double& x, const double& y) {
+  return x * x - y * x + y * y;
+}
+
+int main(int argc, char** argv) {
+  google::InitGoogleLogging(argv[0]);
+  // Problem sizes
+  const int kGridRowsHalf = 9;
+  const int kGridColsHalf = 11;
+  const int kGridRows = 2 * kGridRowsHalf + 1;
+  const int kGridCols = 2 * kGridColsHalf + 1;
+  const int kPoints = 4;
+
+  const Eigen::Vector2d shift(1.234, 2.345);
+  const std::array<Eigen::Vector2d, kPoints> points = {
+      Eigen::Vector2d{-2., -3.},
+      Eigen::Vector2d{-2., 3.},
+      Eigen::Vector2d{2., 3.},
+      Eigen::Vector2d{2., -3.}};
+
+  // Data is a row-major array of kGridRows x kGridCols values of function
+  // f(x, y) on the grid, with x in {-kGridColsHalf, ..., +kGridColsHalf},
+  // and y in {-kGridRowsHalf, ..., +kGridRowsHalf}
+  double data[kGridRows * kGridCols];
+  for (int i = 0; i < kGridRows; ++i) {
+    for (int j = 0; j < kGridCols; ++j) {
+      // Using row-major order
+      int index = i * kGridCols + j;
+      double y = i - kGridRowsHalf;
+      double x = j - kGridColsHalf;
+
+      data[index] = f(x, y);
+    }
+  }
+  const Grid grid(data,
+                  -kGridRowsHalf,
+                  kGridRowsHalf + 1,
+                  -kGridColsHalf,
+                  kGridColsHalf + 1);
+  const Interpolator interpolator(grid);
+
+  Eigen::Vector2d shift_estimate(3.1415, 1.337);
+
+  ceres::Problem problem;
+  problem.AddParameterBlock(shift_estimate.data(), 2);
+
+  for (const auto& p : points) {
+    const Eigen::Vector2d shifted = p + shift;
+
+    const double v = f(shifted.x(), shifted.y());
+    problem.AddResidualBlock(AnalyticBiCubicCost::Create(interpolator, p, v),
+                             nullptr,
+                             shift_estimate.data());
+  }
+
+  ceres::Solver::Options options;
+  options.minimizer_progress_to_stdout = true;
+
+  ceres::Solver::Summary summary;
+  ceres::Solve(options, &problem, &summary);
+  std::cout << summary.BriefReport() << '\n';
+
+  std::cout << "Bicubic interpolation with analytic derivatives:\n";
+  std::cout << "Estimated shift: " << shift_estimate.transpose()
+            << ", ground-truth: " << shift.transpose()
+            << " (error: " << (shift_estimate - shift).transpose() << ")"
+            << std::endl;
+
+  CHECK_LT((shift_estimate - shift).norm(), 1e-9);
+  return 0;
+}
diff --git a/examples/bundle_adjuster.cc b/examples/bundle_adjuster.cc
index e7b154e..582ae2e 100644
--- a/examples/bundle_adjuster.cc
+++ b/examples/bundle_adjuster.cc
@@ -1,5 +1,5 @@
 // Ceres Solver - A fast non-linear least squares minimizer
-// Copyright 2015 Google Inc. All rights reserved.
+// Copyright 2023 Google Inc. All rights reserved.
 // http://ceres-solver.org/
 //
 // Redistribution and use in source and binary forms, with or without
@@ -55,7 +55,9 @@
 #include <cmath>
 #include <cstdio>
 #include <cstdlib>
+#include <memory>
 #include <string>
+#include <thread>
 #include <vector>
 
 #include "bal_problem.h"
@@ -80,34 +82,45 @@
               "automatic, cameras, points, cameras,points, points,cameras");
 
 DEFINE_string(linear_solver, "sparse_schur", "Options are: "
-              "sparse_schur, dense_schur, iterative_schur, sparse_normal_cholesky, "
-              "dense_qr, dense_normal_cholesky and cgnr.");
+              "sparse_schur, dense_schur, iterative_schur, "
+              "sparse_normal_cholesky, dense_qr, dense_normal_cholesky, "
+              "and cgnr.");
 DEFINE_bool(explicit_schur_complement, false, "If using ITERATIVE_SCHUR "
             "then explicitly compute the Schur complement.");
 DEFINE_string(preconditioner, "jacobi", "Options are: "
-              "identity, jacobi, schur_jacobi, cluster_jacobi, "
+              "identity, jacobi, schur_jacobi, schur_power_series_expansion, cluster_jacobi, "
               "cluster_tridiagonal.");
 DEFINE_string(visibility_clustering, "canonical_views",
               "single_linkage, canonical_views");
+DEFINE_bool(use_spse_initialization, false,
+            "Use power series expansion to initialize the solution in ITERATIVE_SCHUR linear solver.");
 
 DEFINE_string(sparse_linear_algebra_library, "suite_sparse",
-              "Options are: suite_sparse and cx_sparse.");
+              "Options are: suite_sparse, accelerate_sparse, eigen_sparse, and "
+              "cuda_sparse.");
 DEFINE_string(dense_linear_algebra_library, "eigen",
-              "Options are: eigen and lapack.");
-DEFINE_string(ordering, "automatic", "Options are: automatic, user.");
+              "Options are: eigen, lapack, and cuda");
+DEFINE_string(ordering_type, "amd", "Options are: amd, nesdis");
+DEFINE_string(linear_solver_ordering, "user",
+              "Options are: automatic and user");
 
 DEFINE_bool(use_quaternions, false, "If true, uses quaternions to represent "
             "rotations. If false, angle axis is used.");
-DEFINE_bool(use_local_parameterization, false, "For quaternions, use a local "
-            "parameterization.");
+DEFINE_bool(use_manifolds, false, "For quaternions, use a manifold.");
 DEFINE_bool(robustify, false, "Use a robust loss function.");
 
 DEFINE_double(eta, 1e-2, "Default value for eta. Eta determines the "
               "accuracy of each linear solve of the truncated newton step. "
               "Changing this parameter can affect solve performance.");
 
-DEFINE_int32(num_threads, 1, "Number of threads.");
+DEFINE_int32(num_threads, -1, "Number of threads. -1 = std::thread::hardware_concurrency.");
 DEFINE_int32(num_iterations, 5, "Number of iterations.");
+DEFINE_int32(max_linear_solver_iterations, 500, "Maximum number of iterations"
+            " for solution of linear system.");
+DEFINE_double(spse_tolerance, 0.1,
+             "Tolerance to reach during the iterations of power series expansion initialization or preconditioning.");
+DEFINE_int32(max_num_spse_iterations, 5,
+             "Maximum number of iterations for power series expansion initialization or preconditioning.");
 DEFINE_double(max_solver_time, 1e32, "Maximum solve time in seconds.");
 DEFINE_bool(nonmonotonic_steps, false, "Trust region algorithm can use"
             " nonmonotic steps.");
@@ -120,7 +133,7 @@
               "perturbation.");
 DEFINE_int32(random_seed, 38401, "Random seed used to set the state "
              "of the pseudo random number generator used to generate "
-             "the pertubations.");
+             "the perturbations.");
 DEFINE_bool(line_search, false, "Use a line search instead of trust region "
             "algorithm.");
 DEFINE_bool(mixed_precision_solves, false, "Use mixed precision solves.");
@@ -128,29 +141,41 @@
 DEFINE_string(initial_ply, "", "Export the BAL file data as a PLY file.");
 DEFINE_string(final_ply, "", "Export the refined BAL file data as a PLY "
               "file.");
-
 // clang-format on
 
-namespace ceres {
-namespace examples {
+namespace ceres::examples {
 namespace {
 
 void SetLinearSolver(Solver::Options* options) {
-  CHECK(StringToLinearSolverType(FLAGS_linear_solver,
+  CHECK(StringToLinearSolverType(CERES_GET_FLAG(FLAGS_linear_solver),
                                  &options->linear_solver_type));
-  CHECK(StringToPreconditionerType(FLAGS_preconditioner,
+  CHECK(StringToPreconditionerType(CERES_GET_FLAG(FLAGS_preconditioner),
                                    &options->preconditioner_type));
-  CHECK(StringToVisibilityClusteringType(FLAGS_visibility_clustering,
-                                         &options->visibility_clustering_type));
+  CHECK(StringToVisibilityClusteringType(
+      CERES_GET_FLAG(FLAGS_visibility_clustering),
+      &options->visibility_clustering_type));
   CHECK(StringToSparseLinearAlgebraLibraryType(
-      FLAGS_sparse_linear_algebra_library,
+      CERES_GET_FLAG(FLAGS_sparse_linear_algebra_library),
       &options->sparse_linear_algebra_library_type));
   CHECK(StringToDenseLinearAlgebraLibraryType(
-      FLAGS_dense_linear_algebra_library,
+      CERES_GET_FLAG(FLAGS_dense_linear_algebra_library),
       &options->dense_linear_algebra_library_type));
-  options->use_explicit_schur_complement = FLAGS_explicit_schur_complement;
-  options->use_mixed_precision_solves = FLAGS_mixed_precision_solves;
-  options->max_num_refinement_iterations = FLAGS_max_num_refinement_iterations;
+  CHECK(
+      StringToLinearSolverOrderingType(CERES_GET_FLAG(FLAGS_ordering_type),
+                                       &options->linear_solver_ordering_type));
+  options->use_explicit_schur_complement =
+      CERES_GET_FLAG(FLAGS_explicit_schur_complement);
+  options->use_mixed_precision_solves =
+      CERES_GET_FLAG(FLAGS_mixed_precision_solves);
+  options->max_num_refinement_iterations =
+      CERES_GET_FLAG(FLAGS_max_num_refinement_iterations);
+  options->max_linear_solver_iterations =
+      CERES_GET_FLAG(FLAGS_max_linear_solver_iterations);
+  options->use_spse_initialization =
+      CERES_GET_FLAG(FLAGS_use_spse_initialization);
+  options->spse_tolerance = CERES_GET_FLAG(FLAGS_spse_tolerance);
+  options->max_num_spse_iterations =
+      CERES_GET_FLAG(FLAGS_max_num_spse_iterations);
 }
 
 void SetOrdering(BALProblem* bal_problem, Solver::Options* options) {
@@ -163,23 +188,27 @@
   double* cameras = bal_problem->mutable_cameras();
 
   if (options->use_inner_iterations) {
-    if (FLAGS_blocks_for_inner_iterations == "cameras") {
+    if (CERES_GET_FLAG(FLAGS_blocks_for_inner_iterations) == "cameras") {
       LOG(INFO) << "Camera blocks for inner iterations";
-      options->inner_iteration_ordering.reset(new ParameterBlockOrdering);
+      options->inner_iteration_ordering =
+          std::make_shared<ParameterBlockOrdering>();
       for (int i = 0; i < num_cameras; ++i) {
         options->inner_iteration_ordering->AddElementToGroup(
             cameras + camera_block_size * i, 0);
       }
-    } else if (FLAGS_blocks_for_inner_iterations == "points") {
+    } else if (CERES_GET_FLAG(FLAGS_blocks_for_inner_iterations) == "points") {
       LOG(INFO) << "Point blocks for inner iterations";
-      options->inner_iteration_ordering.reset(new ParameterBlockOrdering);
+      options->inner_iteration_ordering =
+          std::make_shared<ParameterBlockOrdering>();
       for (int i = 0; i < num_points; ++i) {
         options->inner_iteration_ordering->AddElementToGroup(
             points + point_block_size * i, 0);
       }
-    } else if (FLAGS_blocks_for_inner_iterations == "cameras,points") {
+    } else if (CERES_GET_FLAG(FLAGS_blocks_for_inner_iterations) ==
+               "cameras,points") {
       LOG(INFO) << "Camera followed by point blocks for inner iterations";
-      options->inner_iteration_ordering.reset(new ParameterBlockOrdering);
+      options->inner_iteration_ordering =
+          std::make_shared<ParameterBlockOrdering>();
       for (int i = 0; i < num_cameras; ++i) {
         options->inner_iteration_ordering->AddElementToGroup(
             cameras + camera_block_size * i, 0);
@@ -188,9 +217,11 @@
         options->inner_iteration_ordering->AddElementToGroup(
             points + point_block_size * i, 1);
       }
-    } else if (FLAGS_blocks_for_inner_iterations == "points,cameras") {
+    } else if (CERES_GET_FLAG(FLAGS_blocks_for_inner_iterations) ==
+               "points,cameras") {
       LOG(INFO) << "Point followed by camera blocks for inner iterations";
-      options->inner_iteration_ordering.reset(new ParameterBlockOrdering);
+      options->inner_iteration_ordering =
+          std::make_shared<ParameterBlockOrdering>();
       for (int i = 0; i < num_cameras; ++i) {
         options->inner_iteration_ordering->AddElementToGroup(
             cameras + camera_block_size * i, 1);
@@ -199,11 +230,12 @@
         options->inner_iteration_ordering->AddElementToGroup(
             points + point_block_size * i, 0);
       }
-    } else if (FLAGS_blocks_for_inner_iterations == "automatic") {
+    } else if (CERES_GET_FLAG(FLAGS_blocks_for_inner_iterations) ==
+               "automatic") {
       LOG(INFO) << "Choosing automatic blocks for inner iterations";
     } else {
       LOG(FATAL) << "Unknown block type for inner iterations: "
-                 << FLAGS_blocks_for_inner_iterations;
+                 << CERES_GET_FLAG(FLAGS_blocks_for_inner_iterations);
     }
   }
 
@@ -213,46 +245,54 @@
   // ITERATIVE_SCHUR solvers make use of this specialized
   // structure.
   //
-  // This can either be done by specifying Options::ordering_type =
-  // ceres::SCHUR, in which case Ceres will automatically determine
-  // the right ParameterBlock ordering, or by manually specifying a
-  // suitable ordering vector and defining
-  // Options::num_eliminate_blocks.
-  if (FLAGS_ordering == "automatic") {
-    return;
+  // This can either be done by specifying a
+  // Options::linear_solver_ordering or having Ceres figure it out
+  // automatically using a greedy maximum independent set algorithm.
+  if (CERES_GET_FLAG(FLAGS_linear_solver_ordering) == "user") {
+    auto* ordering = new ceres::ParameterBlockOrdering;
+
+    // The points come before the cameras.
+    for (int i = 0; i < num_points; ++i) {
+      ordering->AddElementToGroup(points + point_block_size * i, 0);
+    }
+
+    for (int i = 0; i < num_cameras; ++i) {
+      // When using axis-angle, there is a single parameter block for
+      // the entire camera.
+      ordering->AddElementToGroup(cameras + camera_block_size * i, 1);
+    }
+
+    options->linear_solver_ordering.reset(ordering);
   }
-
-  ceres::ParameterBlockOrdering* ordering = new ceres::ParameterBlockOrdering;
-
-  // The points come before the cameras.
-  for (int i = 0; i < num_points; ++i) {
-    ordering->AddElementToGroup(points + point_block_size * i, 0);
-  }
-
-  for (int i = 0; i < num_cameras; ++i) {
-    // When using axis-angle, there is a single parameter block for
-    // the entire camera.
-    ordering->AddElementToGroup(cameras + camera_block_size * i, 1);
-  }
-
-  options->linear_solver_ordering.reset(ordering);
 }
 
 void SetMinimizerOptions(Solver::Options* options) {
-  options->max_num_iterations = FLAGS_num_iterations;
+  options->max_num_iterations = CERES_GET_FLAG(FLAGS_num_iterations);
   options->minimizer_progress_to_stdout = true;
-  options->num_threads = FLAGS_num_threads;
-  options->eta = FLAGS_eta;
-  options->max_solver_time_in_seconds = FLAGS_max_solver_time;
-  options->use_nonmonotonic_steps = FLAGS_nonmonotonic_steps;
-  if (FLAGS_line_search) {
+  if (CERES_GET_FLAG(FLAGS_num_threads) == -1) {
+    const int num_available_threads =
+        static_cast<int>(std::thread::hardware_concurrency());
+    if (num_available_threads > 0) {
+      options->num_threads = num_available_threads;
+    }
+  } else {
+    options->num_threads = CERES_GET_FLAG(FLAGS_num_threads);
+  }
+  CHECK_GE(options->num_threads, 1);
+
+  options->eta = CERES_GET_FLAG(FLAGS_eta);
+  options->max_solver_time_in_seconds = CERES_GET_FLAG(FLAGS_max_solver_time);
+  options->use_nonmonotonic_steps = CERES_GET_FLAG(FLAGS_nonmonotonic_steps);
+  if (CERES_GET_FLAG(FLAGS_line_search)) {
     options->minimizer_type = ceres::LINE_SEARCH;
   }
 
-  CHECK(StringToTrustRegionStrategyType(FLAGS_trust_region_strategy,
-                                        &options->trust_region_strategy_type));
-  CHECK(StringToDoglegType(FLAGS_dogleg, &options->dogleg_type));
-  options->use_inner_iterations = FLAGS_inner_iterations;
+  CHECK(StringToTrustRegionStrategyType(
+      CERES_GET_FLAG(FLAGS_trust_region_strategy),
+      &options->trust_region_strategy_type));
+  CHECK(
+      StringToDoglegType(CERES_GET_FLAG(FLAGS_dogleg), &options->dogleg_type));
+  options->use_inner_iterations = CERES_GET_FLAG(FLAGS_inner_iterations);
 }
 
 void SetSolverOptionsFromFlags(BALProblem* bal_problem,
@@ -276,16 +316,17 @@
     CostFunction* cost_function;
     // Each Residual block takes a point and a camera as input and
     // outputs a 2 dimensional residual.
-    cost_function = (FLAGS_use_quaternions)
+    cost_function = (CERES_GET_FLAG(FLAGS_use_quaternions))
                         ? SnavelyReprojectionErrorWithQuaternions::Create(
                               observations[2 * i + 0], observations[2 * i + 1])
                         : SnavelyReprojectionError::Create(
                               observations[2 * i + 0], observations[2 * i + 1]);
 
     // If enabled use Huber's loss function.
-    LossFunction* loss_function = FLAGS_robustify ? new HuberLoss(1.0) : NULL;
+    LossFunction* loss_function =
+        CERES_GET_FLAG(FLAGS_robustify) ? new HuberLoss(1.0) : nullptr;
 
-    // Each observation correponds to a pair of a camera and a point
+    // Each observation corresponds to a pair of a camera and a point
     // which are identified by camera_index()[i] and point_index()[i]
     // respectively.
     double* camera =
@@ -294,60 +335,60 @@
     problem->AddResidualBlock(cost_function, loss_function, camera, point);
   }
 
-  if (FLAGS_use_quaternions && FLAGS_use_local_parameterization) {
-    LocalParameterization* camera_parameterization =
-        new ProductParameterization(new QuaternionParameterization(),
-                                    new IdentityParameterization(6));
+  if (CERES_GET_FLAG(FLAGS_use_quaternions) &&
+      CERES_GET_FLAG(FLAGS_use_manifolds)) {
+    Manifold* camera_manifold =
+        new ProductManifold<QuaternionManifold, EuclideanManifold<6>>{};
     for (int i = 0; i < bal_problem->num_cameras(); ++i) {
-      problem->SetParameterization(cameras + camera_block_size * i,
-                                   camera_parameterization);
+      problem->SetManifold(cameras + camera_block_size * i, camera_manifold);
     }
   }
 }
 
 void SolveProblem(const char* filename) {
-  BALProblem bal_problem(filename, FLAGS_use_quaternions);
+  BALProblem bal_problem(filename, CERES_GET_FLAG(FLAGS_use_quaternions));
 
-  if (!FLAGS_initial_ply.empty()) {
-    bal_problem.WriteToPLYFile(FLAGS_initial_ply);
+  if (!CERES_GET_FLAG(FLAGS_initial_ply).empty()) {
+    bal_problem.WriteToPLYFile(CERES_GET_FLAG(FLAGS_initial_ply));
   }
 
   Problem problem;
 
-  srand(FLAGS_random_seed);
+  srand(CERES_GET_FLAG(FLAGS_random_seed));
   bal_problem.Normalize();
-  bal_problem.Perturb(
-      FLAGS_rotation_sigma, FLAGS_translation_sigma, FLAGS_point_sigma);
+  bal_problem.Perturb(CERES_GET_FLAG(FLAGS_rotation_sigma),
+                      CERES_GET_FLAG(FLAGS_translation_sigma),
+                      CERES_GET_FLAG(FLAGS_point_sigma));
 
   BuildProblem(&bal_problem, &problem);
   Solver::Options options;
   SetSolverOptionsFromFlags(&bal_problem, &options);
   options.gradient_tolerance = 1e-16;
   options.function_tolerance = 1e-16;
+  options.parameter_tolerance = 1e-16;
   Solver::Summary summary;
   Solve(options, &problem, &summary);
   std::cout << summary.FullReport() << "\n";
 
-  if (!FLAGS_final_ply.empty()) {
-    bal_problem.WriteToPLYFile(FLAGS_final_ply);
+  if (!CERES_GET_FLAG(FLAGS_final_ply).empty()) {
+    bal_problem.WriteToPLYFile(CERES_GET_FLAG(FLAGS_final_ply));
   }
 }
 
 }  // namespace
-}  // namespace examples
-}  // namespace ceres
+}  // namespace ceres::examples
 
 int main(int argc, char** argv) {
   GFLAGS_NAMESPACE::ParseCommandLineFlags(&argc, &argv, true);
   google::InitGoogleLogging(argv[0]);
-  if (FLAGS_input.empty()) {
+  if (CERES_GET_FLAG(FLAGS_input).empty()) {
     LOG(ERROR) << "Usage: bundle_adjuster --input=bal_problem";
     return 1;
   }
 
-  CHECK(FLAGS_use_quaternions || !FLAGS_use_local_parameterization)
-      << "--use_local_parameterization can only be used with "
-      << "--use_quaternions.";
-  ceres::examples::SolveProblem(FLAGS_input.c_str());
+  CHECK(CERES_GET_FLAG(FLAGS_use_quaternions) ||
+        !CERES_GET_FLAG(FLAGS_use_manifolds))
+      << "--use_manifolds can only be used with --use_quaternions.";
+  ceres::examples::SolveProblem(CERES_GET_FLAG(FLAGS_input).c_str());
   return 0;
 }
diff --git a/examples/circle_fit.cc b/examples/circle_fit.cc
index c542475..fd848d9 100644
--- a/examples/circle_fit.cc
+++ b/examples/circle_fit.cc
@@ -1,5 +1,5 @@
 // Ceres Solver - A fast non-linear least squares minimizer
-// Copyright 2015 Google Inc. All rights reserved.
+// Copyright 2023 Google Inc. All rights reserved.
 // http://ceres-solver.org/
 //
 // Redistribution and use in source and binary forms, with or without
@@ -57,14 +57,6 @@
 #include "gflags/gflags.h"
 #include "glog/logging.h"
 
-using ceres::AutoDiffCostFunction;
-using ceres::CauchyLoss;
-using ceres::CostFunction;
-using ceres::LossFunction;
-using ceres::Problem;
-using ceres::Solve;
-using ceres::Solver;
-
 DEFINE_double(robust_threshold,
               0.0,
               "Robust loss parameter. Set to 0 for normal squared error (no "
@@ -128,21 +120,21 @@
   // Parameterize r as m^2 so that it can't be negative.
   double m = sqrt(r);
 
-  Problem problem;
+  ceres::Problem problem;
 
   // Configure the loss function.
-  LossFunction* loss = NULL;
-  if (FLAGS_robust_threshold) {
-    loss = new CauchyLoss(FLAGS_robust_threshold);
+  ceres::LossFunction* loss = nullptr;
+  if (CERES_GET_FLAG(FLAGS_robust_threshold)) {
+    loss = new ceres::CauchyLoss(CERES_GET_FLAG(FLAGS_robust_threshold));
   }
 
   // Add the residuals.
   double xx, yy;
   int num_points = 0;
   while (scanf("%lf %lf\n", &xx, &yy) == 2) {
-    CostFunction* cost =
-        new AutoDiffCostFunction<DistanceFromCircleCost, 1, 1, 1, 1>(
-            new DistanceFromCircleCost(xx, yy));
+    ceres::CostFunction* cost =
+        new ceres::AutoDiffCostFunction<DistanceFromCircleCost, 1, 1, 1, 1>(xx,
+                                                                            yy);
     problem.AddResidualBlock(cost, loss, &x, &y, &m);
     num_points++;
   }
@@ -150,11 +142,11 @@
   std::cout << "Got " << num_points << " points.\n";
 
   // Build and solve the problem.
-  Solver::Options options;
+  ceres::Solver::Options options;
   options.max_num_iterations = 500;
   options.linear_solver_type = ceres::DENSE_QR;
-  Solver::Summary summary;
-  Solve(options, &problem, &summary);
+  ceres::Solver::Summary summary;
+  ceres::Solve(options, &problem, &summary);
 
   // Recover r from m.
   r = m * m;
diff --git a/examples/curve_fitting.cc b/examples/curve_fitting.cc
index fc7ff94..105402e 100644
--- a/examples/curve_fitting.cc
+++ b/examples/curve_fitting.cc
@@ -1,5 +1,5 @@
 // Ceres Solver - A fast non-linear least squares minimizer
-// Copyright 2015 Google Inc. All rights reserved.
+// Copyright 2023 Google Inc. All rights reserved.
 // http://ceres-solver.org/
 //
 // Redistribution and use in source and binary forms, with or without
@@ -27,16 +27,13 @@
 // POSSIBILITY OF SUCH DAMAGE.
 //
 // Author: sameeragarwal@google.com (Sameer Agarwal)
+//
+// This example fits the curve f(x;m,c) = e^(m * x + c) to data, minimizing the
+// sum squared loss.
 
 #include "ceres/ceres.h"
 #include "glog/logging.h"
 
-using ceres::AutoDiffCostFunction;
-using ceres::CostFunction;
-using ceres::Problem;
-using ceres::Solve;
-using ceres::Solver;
-
 // Data generated using the following octave code.
 //   randn('seed', 23497);
 //   m = 0.3;
@@ -137,28 +134,30 @@
 int main(int argc, char** argv) {
   google::InitGoogleLogging(argv[0]);
 
-  double m = 0.0;
-  double c = 0.0;
+  const double initial_m = 0.0;
+  const double initial_c = 0.0;
+  double m = initial_m;
+  double c = initial_c;
 
-  Problem problem;
+  ceres::Problem problem;
   for (int i = 0; i < kNumObservations; ++i) {
     problem.AddResidualBlock(
-        new AutoDiffCostFunction<ExponentialResidual, 1, 1, 1>(
-            new ExponentialResidual(data[2 * i], data[2 * i + 1])),
-        NULL,
+        new ceres::AutoDiffCostFunction<ExponentialResidual, 1, 1, 1>(
+            data[2 * i], data[2 * i + 1]),
+        nullptr,
         &m,
         &c);
   }
 
-  Solver::Options options;
+  ceres::Solver::Options options;
   options.max_num_iterations = 25;
   options.linear_solver_type = ceres::DENSE_QR;
   options.minimizer_progress_to_stdout = true;
 
-  Solver::Summary summary;
-  Solve(options, &problem, &summary);
+  ceres::Solver::Summary summary;
+  ceres::Solve(options, &problem, &summary);
   std::cout << summary.BriefReport() << "\n";
-  std::cout << "Initial m: " << 0.0 << " c: " << 0.0 << "\n";
+  std::cout << "Initial m: " << initial_m << " c: " << initial_c << "\n";
   std::cout << "Final   m: " << m << " c: " << c << "\n";
   return 0;
 }
diff --git a/examples/denoising.cc b/examples/denoising.cc
index 61ea2c6..dc13d19 100644
--- a/examples/denoising.cc
+++ b/examples/denoising.cc
@@ -1,5 +1,5 @@
 // Ceres Solver - A fast non-linear least squares minimizer
-// Copyright 2015 Google Inc. All rights reserved.
+// Copyright 2023 Google Inc. All rights reserved.
 // http://ceres-solver.org/
 //
 // Redistribution and use in source and binary forms, with or without
@@ -33,7 +33,7 @@
 // Note that for good denoising results the weighting between the data term
 // and the Fields of Experts term needs to be adjusted. This is discussed
 // in [1]. This program assumes Gaussian noise. The noise model can be changed
-// by substituing another function for QuadraticCostFunction.
+// by substituting another function for QuadraticCostFunction.
 //
 // [1] S. Roth and M.J. Black. "Fields of Experts." International Journal of
 //     Computer Vision, 82(2):205--229, 2009.
@@ -102,8 +102,7 @@
               "The fraction of residual blocks to use for the"
               " subset preconditioner.");
 
-namespace ceres {
-namespace examples {
+namespace ceres::examples {
 namespace {
 
 // This cost function is used to build the data term.
@@ -113,12 +112,12 @@
 class QuadraticCostFunction : public ceres::SizedCostFunction<1, 1> {
  public:
   QuadraticCostFunction(double a, double b) : sqrta_(std::sqrt(a)), b_(b) {}
-  virtual bool Evaluate(double const* const* parameters,
-                        double* residuals,
-                        double** jacobians) const {
+  bool Evaluate(double const* const* parameters,
+                double* residuals,
+                double** jacobians) const override {
     const double x = parameters[0][0];
     residuals[0] = sqrta_ * (x - b_);
-    if (jacobians != NULL && jacobians[0] != NULL) {
+    if (jacobians != nullptr && jacobians[0] != nullptr) {
       jacobians[0][0] = sqrta_;
     }
     return true;
@@ -134,13 +133,14 @@
                    Problem* problem,
                    PGMImage<double>* solution) {
   // Create the data term
-  CHECK_GT(FLAGS_sigma, 0.0);
-  const double coefficient = 1 / (2.0 * FLAGS_sigma * FLAGS_sigma);
+  CHECK_GT(CERES_GET_FLAG(FLAGS_sigma), 0.0);
+  const double coefficient =
+      1 / (2.0 * CERES_GET_FLAG(FLAGS_sigma) * CERES_GET_FLAG(FLAGS_sigma));
   for (int index = 0; index < image.NumPixels(); ++index) {
     ceres::CostFunction* cost_function = new QuadraticCostFunction(
         coefficient, image.PixelFromLinearIndex(index));
     problem->AddResidualBlock(
-        cost_function, NULL, solution->MutablePixelFromLinearIndex(index));
+        cost_function, nullptr, solution->MutablePixelFromLinearIndex(index));
   }
 
   // Create Ceres cost and loss functions for regularization. One is needed for
@@ -175,31 +175,35 @@
 }
 
 void SetLinearSolver(Solver::Options* options) {
-  CHECK(StringToLinearSolverType(FLAGS_linear_solver,
+  CHECK(StringToLinearSolverType(CERES_GET_FLAG(FLAGS_linear_solver),
                                  &options->linear_solver_type));
-  CHECK(StringToPreconditionerType(FLAGS_preconditioner,
+  CHECK(StringToPreconditionerType(CERES_GET_FLAG(FLAGS_preconditioner),
                                    &options->preconditioner_type));
   CHECK(StringToSparseLinearAlgebraLibraryType(
-      FLAGS_sparse_linear_algebra_library,
+      CERES_GET_FLAG(FLAGS_sparse_linear_algebra_library),
       &options->sparse_linear_algebra_library_type));
-  options->use_mixed_precision_solves = FLAGS_mixed_precision_solves;
-  options->max_num_refinement_iterations = FLAGS_max_num_refinement_iterations;
+  options->use_mixed_precision_solves =
+      CERES_GET_FLAG(FLAGS_mixed_precision_solves);
+  options->max_num_refinement_iterations =
+      CERES_GET_FLAG(FLAGS_max_num_refinement_iterations);
 }
 
 void SetMinimizerOptions(Solver::Options* options) {
-  options->max_num_iterations = FLAGS_num_iterations;
+  options->max_num_iterations = CERES_GET_FLAG(FLAGS_num_iterations);
   options->minimizer_progress_to_stdout = true;
-  options->num_threads = FLAGS_num_threads;
-  options->eta = FLAGS_eta;
-  options->use_nonmonotonic_steps = FLAGS_nonmonotonic_steps;
-  if (FLAGS_line_search) {
+  options->num_threads = CERES_GET_FLAG(FLAGS_num_threads);
+  options->eta = CERES_GET_FLAG(FLAGS_eta);
+  options->use_nonmonotonic_steps = CERES_GET_FLAG(FLAGS_nonmonotonic_steps);
+  if (CERES_GET_FLAG(FLAGS_line_search)) {
     options->minimizer_type = ceres::LINE_SEARCH;
   }
 
-  CHECK(StringToTrustRegionStrategyType(FLAGS_trust_region_strategy,
-                                        &options->trust_region_strategy_type));
-  CHECK(StringToDoglegType(FLAGS_dogleg, &options->dogleg_type));
-  options->use_inner_iterations = FLAGS_inner_iterations;
+  CHECK(StringToTrustRegionStrategyType(
+      CERES_GET_FLAG(FLAGS_trust_region_strategy),
+      &options->trust_region_strategy_type));
+  CHECK(
+      StringToDoglegType(CERES_GET_FLAG(FLAGS_dogleg), &options->dogleg_type));
+  options->use_inner_iterations = CERES_GET_FLAG(FLAGS_inner_iterations);
 }
 
 // Solves the FoE problem using Ceres and post-processes it to make sure the
@@ -226,7 +230,7 @@
     std::default_random_engine engine;
     std::uniform_real_distribution<> distribution(0, 1);  // rage 0 - 1
     for (auto residual_block : residual_blocks) {
-      if (distribution(engine) <= FLAGS_subset_fraction) {
+      if (distribution(engine) <= CERES_GET_FLAG(FLAGS_subset_fraction)) {
         options.residual_blocks_for_subset_preconditioner.insert(
             residual_block);
       }
@@ -247,20 +251,19 @@
 }
 
 }  // namespace
-}  // namespace examples
-}  // namespace ceres
+}  // namespace ceres::examples
 
 int main(int argc, char** argv) {
   using namespace ceres::examples;
   GFLAGS_NAMESPACE::ParseCommandLineFlags(&argc, &argv, true);
   google::InitGoogleLogging(argv[0]);
 
-  if (FLAGS_input.empty()) {
+  if (CERES_GET_FLAG(FLAGS_input).empty()) {
     std::cerr << "Please provide an image file name using -input.\n";
     return 1;
   }
 
-  if (FLAGS_foe_file.empty()) {
+  if (CERES_GET_FLAG(FLAGS_foe_file).empty()) {
     std::cerr << "Please provide a Fields of Experts file name using -foe_file."
                  "\n";
     return 1;
@@ -268,15 +271,16 @@
 
   // Load the Fields of Experts filters from file.
   FieldsOfExperts foe;
-  if (!foe.LoadFromFile(FLAGS_foe_file)) {
-    std::cerr << "Loading \"" << FLAGS_foe_file << "\" failed.\n";
+  if (!foe.LoadFromFile(CERES_GET_FLAG(FLAGS_foe_file))) {
+    std::cerr << "Loading \"" << CERES_GET_FLAG(FLAGS_foe_file)
+              << "\" failed.\n";
     return 2;
   }
 
   // Read the images
-  PGMImage<double> image(FLAGS_input);
+  PGMImage<double> image(CERES_GET_FLAG(FLAGS_input));
   if (image.width() == 0) {
-    std::cerr << "Reading \"" << FLAGS_input << "\" failed.\n";
+    std::cerr << "Reading \"" << CERES_GET_FLAG(FLAGS_input) << "\" failed.\n";
     return 3;
   }
   PGMImage<double> solution(image.width(), image.height());
@@ -287,9 +291,9 @@
 
   SolveProblem(&problem, &solution);
 
-  if (!FLAGS_output.empty()) {
-    CHECK(solution.WriteToFile(FLAGS_output))
-        << "Writing \"" << FLAGS_output << "\" failed.";
+  if (!CERES_GET_FLAG(FLAGS_output).empty()) {
+    CHECK(solution.WriteToFile(CERES_GET_FLAG(FLAGS_output)))
+        << "Writing \"" << CERES_GET_FLAG(FLAGS_output) << "\" failed.";
   }
 
   return 0;
diff --git a/examples/ellipse_approximation.cc b/examples/ellipse_approximation.cc
index 74782f4..6fa8f1c 100644
--- a/examples/ellipse_approximation.cc
+++ b/examples/ellipse_approximation.cc
@@ -1,5 +1,5 @@
 // Ceres Solver - A fast non-linear least squares minimizer
-// Copyright 2015 Google Inc. All rights reserved.
+// Copyright 2023 Google Inc. All rights reserved.
 // http://ceres-solver.org/
 //
 // Redistribution and use in source and binary forms, with or without
@@ -36,6 +36,7 @@
 // dense but dynamically sparse.
 
 #include <cmath>
+#include <utility>
 #include <vector>
 
 #include "ceres/ceres.h"
@@ -275,8 +276,8 @@
 class PointToLineSegmentContourCostFunction : public ceres::CostFunction {
  public:
   PointToLineSegmentContourCostFunction(const int num_segments,
-                                        const Eigen::Vector2d& y)
-      : num_segments_(num_segments), y_(y) {
+                                        Eigen::Vector2d y)
+      : num_segments_(num_segments), y_(std::move(y)) {
     // The first parameter is the preimage position.
     mutable_parameter_block_sizes()->push_back(1);
     // The next parameters are the control points for the line segment contour.
@@ -286,9 +287,9 @@
     set_num_residuals(2);
   }
 
-  virtual bool Evaluate(const double* const* x,
-                        double* residuals,
-                        double** jacobians) const {
+  bool Evaluate(const double* const* x,
+                double* residuals,
+                double** jacobians) const override {
     // Convert the preimage position `t` into a segment index `i0` and the
     // line segment interpolation parameter `u`. `i1` is the index of the next
     // control point.
@@ -302,16 +303,16 @@
     residuals[0] = y_[0] - ((1.0 - u) * x[1 + i0][0] + u * x[1 + i1][0]);
     residuals[1] = y_[1] - ((1.0 - u) * x[1 + i0][1] + u * x[1 + i1][1]);
 
-    if (jacobians == NULL) {
+    if (jacobians == nullptr) {
       return true;
     }
 
-    if (jacobians[0] != NULL) {
+    if (jacobians[0] != nullptr) {
       jacobians[0][0] = x[1 + i0][0] - x[1 + i1][0];
       jacobians[0][1] = x[1 + i0][1] - x[1 + i1][1];
     }
     for (int i = 0; i < num_segments_; ++i) {
-      if (jacobians[i + 1] != NULL) {
+      if (jacobians[i + 1] != nullptr) {
         ceres::MatrixRef(jacobians[i + 1], 2, 2).setZero();
         if (i == i0) {
           jacobians[i + 1][0] = -(1.0 - u);
@@ -353,7 +354,7 @@
 
   static ceres::CostFunction* Create(const double sqrt_weight) {
     return new ceres::AutoDiffCostFunction<EuclideanDistanceFunctor, 2, 2, 2>(
-        new EuclideanDistanceFunctor(sqrt_weight));
+        sqrt_weight);
   }
 
  private:
@@ -385,8 +386,8 @@
 
   // Eigen::MatrixXd is column major so we define our own MatrixXd which is
   // row major. Eigen::VectorXd can be used directly.
-  typedef Eigen::Matrix<double, Eigen::Dynamic, Eigen::Dynamic, Eigen::RowMajor>
-      MatrixXd;
+  using MatrixXd =
+      Eigen::Matrix<double, Eigen::Dynamic, Eigen::Dynamic, Eigen::RowMajor>;
   using Eigen::VectorXd;
 
   // `X` is the matrix of control points which make up the contour of line
@@ -395,7 +396,7 @@
   //
   // Initialize `X` to points on the unit circle.
   VectorXd w(num_segments + 1);
-  w.setLinSpaced(num_segments + 1, 0.0, 2.0 * M_PI);
+  w.setLinSpaced(num_segments + 1, 0.0, 2.0 * ceres::constants::pi);
   w.conservativeResize(num_segments);
   MatrixXd X(num_segments, 2);
   X.col(0) = w.array().cos();
@@ -404,9 +405,9 @@
   // Each data point has an associated preimage position on the line segment
   // contour. For each data point we initialize the preimage positions to
   // the index of the closest control point.
-  const int num_observations = kY.rows();
+  const int64_t num_observations = kY.rows();
   VectorXd t(num_observations);
-  for (int i = 0; i < num_observations; ++i) {
+  for (int64_t i = 0; i < num_observations; ++i) {
     (X.rowwise() - kY.row(i)).rowwise().squaredNorm().minCoeff(&t[i]);
   }
 
@@ -415,7 +416,7 @@
   // For each data point add a residual which measures its distance to its
   // corresponding position on the line segment contour.
   std::vector<double*> parameter_blocks(1 + num_segments);
-  parameter_blocks[0] = NULL;
+  parameter_blocks[0] = nullptr;
   for (int i = 0; i < num_segments; ++i) {
     parameter_blocks[i + 1] = X.data() + 2 * i;
   }
@@ -423,7 +424,7 @@
     parameter_blocks[0] = &t[i];
     problem.AddResidualBlock(
         PointToLineSegmentContourCostFunction::Create(num_segments, kY.row(i)),
-        NULL,
+        nullptr,
         parameter_blocks);
   }
 
@@ -431,7 +432,7 @@
   for (int i = 0; i < num_segments; ++i) {
     problem.AddResidualBlock(
         EuclideanDistanceFunctor::Create(sqrt(regularization_weight)),
-        NULL,
+        nullptr,
         X.data() + 2 * i,
         X.data() + 2 * ((i + 1) % num_segments));
   }
diff --git a/examples/evaluation_callback_example.cc b/examples/evaluation_callback_example.cc
new file mode 100644
index 0000000..6dbf932
--- /dev/null
+++ b/examples/evaluation_callback_example.cc
@@ -0,0 +1,257 @@
+// Ceres Solver - A fast non-linear least squares minimizer
+// Copyright 2023 Google Inc. All rights reserved.
+// http://ceres-solver.org/
+//
+// Redistribution and use in source and binary forms, with or without
+// modification, are permitted provided that the following conditions are met:
+//
+// * Redistributions of source code must retain the above copyright notice,
+//   this list of conditions and the following disclaimer.
+// * Redistributions in binary form must reproduce the above copyright notice,
+//   this list of conditions and the following disclaimer in the documentation
+//   and/or other materials provided with the distribution.
+// * Neither the name of Google Inc. nor the names of its contributors may be
+//   used to endorse or promote products derived from this software without
+//   specific prior written permission.
+//
+// THIS SOFTWARE IS PROVIDED BY THE COPYRIGHT HOLDERS AND CONTRIBUTORS "AS IS"
+// AND ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED TO, THE
+// IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR PURPOSE
+// ARE DISCLAIMED. IN NO EVENT SHALL THE COPYRIGHT OWNER OR CONTRIBUTORS BE
+// LIABLE FOR ANY DIRECT, INDIRECT, INCIDENTAL, SPECIAL, EXEMPLARY, OR
+// CONSEQUENTIAL DAMAGES (INCLUDING, BUT NOT LIMITED TO, PROCUREMENT OF
+// SUBSTITUTE GOODS OR SERVICES; LOSS OF USE, DATA, OR PROFITS; OR BUSINESS
+// INTERRUPTION) HOWEVER CAUSED AND ON ANY THEORY OF LIABILITY, WHETHER IN
+// CONTRACT, STRICT LIABILITY, OR TORT (INCLUDING NEGLIGENCE OR OTHERWISE)
+// ARISING IN ANY WAY OUT OF THE USE OF THIS SOFTWARE, EVEN IF ADVISED OF THE
+// POSSIBILITY OF SUCH DAMAGE.
+//
+// Author: sameeragarwal@google.com (Sameer Agarwal)
+//
+// This example illustrates the use of the EvaluationCallback, which can be used
+// to perform high performance computation of the residual and Jacobians outside
+// Ceres (in this case using Eigen's vectorized code) and then the CostFunctions
+// just copy these computed residuals and Jacobians appropriately and pass them
+// to Ceres Solver.
+//
+// The results of running this example should be identical to the results
+// obtained by running curve_fitting.cc. The only difference between the two
+// examples is how the residuals and Jacobians are computed.
+//
+// The observant reader will note that both here and curve_fitting.cc instead of
+// creating one ResidualBlock for each observation one can just do one
+// ResidualBlock/CostFunction for the entire problem. The reason for keeping one
+// residual per observation is that it is what is needed if and when we need to
+// introduce a loss function which is what we do in robust_curve_fitting.cc
+
+#include <iostream>
+
+#include "Eigen/Core"
+#include "ceres/ceres.h"
+#include "glog/logging.h"
+
+// Data generated using the following octave code.
+//   randn('seed', 23497);
+//   m = 0.3;
+//   c = 0.1;
+//   x=[0:0.075:5];
+//   y = exp(m * x + c);
+//   noise = randn(size(x)) * 0.2;
+//   y_observed = y + noise;
+//   data = [x', y_observed'];
+
+const int kNumObservations = 67;
+// clang-format off
+const double data[] = {
+  0.000000e+00, 1.133898e+00,
+  7.500000e-02, 1.334902e+00,
+  1.500000e-01, 1.213546e+00,
+  2.250000e-01, 1.252016e+00,
+  3.000000e-01, 1.392265e+00,
+  3.750000e-01, 1.314458e+00,
+  4.500000e-01, 1.472541e+00,
+  5.250000e-01, 1.536218e+00,
+  6.000000e-01, 1.355679e+00,
+  6.750000e-01, 1.463566e+00,
+  7.500000e-01, 1.490201e+00,
+  8.250000e-01, 1.658699e+00,
+  9.000000e-01, 1.067574e+00,
+  9.750000e-01, 1.464629e+00,
+  1.050000e+00, 1.402653e+00,
+  1.125000e+00, 1.713141e+00,
+  1.200000e+00, 1.527021e+00,
+  1.275000e+00, 1.702632e+00,
+  1.350000e+00, 1.423899e+00,
+  1.425000e+00, 1.543078e+00,
+  1.500000e+00, 1.664015e+00,
+  1.575000e+00, 1.732484e+00,
+  1.650000e+00, 1.543296e+00,
+  1.725000e+00, 1.959523e+00,
+  1.800000e+00, 1.685132e+00,
+  1.875000e+00, 1.951791e+00,
+  1.950000e+00, 2.095346e+00,
+  2.025000e+00, 2.361460e+00,
+  2.100000e+00, 2.169119e+00,
+  2.175000e+00, 2.061745e+00,
+  2.250000e+00, 2.178641e+00,
+  2.325000e+00, 2.104346e+00,
+  2.400000e+00, 2.584470e+00,
+  2.475000e+00, 1.914158e+00,
+  2.550000e+00, 2.368375e+00,
+  2.625000e+00, 2.686125e+00,
+  2.700000e+00, 2.712395e+00,
+  2.775000e+00, 2.499511e+00,
+  2.850000e+00, 2.558897e+00,
+  2.925000e+00, 2.309154e+00,
+  3.000000e+00, 2.869503e+00,
+  3.075000e+00, 3.116645e+00,
+  3.150000e+00, 3.094907e+00,
+  3.225000e+00, 2.471759e+00,
+  3.300000e+00, 3.017131e+00,
+  3.375000e+00, 3.232381e+00,
+  3.450000e+00, 2.944596e+00,
+  3.525000e+00, 3.385343e+00,
+  3.600000e+00, 3.199826e+00,
+  3.675000e+00, 3.423039e+00,
+  3.750000e+00, 3.621552e+00,
+  3.825000e+00, 3.559255e+00,
+  3.900000e+00, 3.530713e+00,
+  3.975000e+00, 3.561766e+00,
+  4.050000e+00, 3.544574e+00,
+  4.125000e+00, 3.867945e+00,
+  4.200000e+00, 4.049776e+00,
+  4.275000e+00, 3.885601e+00,
+  4.350000e+00, 4.110505e+00,
+  4.425000e+00, 4.345320e+00,
+  4.500000e+00, 4.161241e+00,
+  4.575000e+00, 4.363407e+00,
+  4.650000e+00, 4.161576e+00,
+  4.725000e+00, 4.619728e+00,
+  4.800000e+00, 4.737410e+00,
+  4.875000e+00, 4.727863e+00,
+  4.950000e+00, 4.669206e+00,
+};
+// clang-format on
+
+// This implementation of the EvaluationCallback interface also stores the
+// residuals and Jacobians that the CostFunction copies their values from.
+class MyEvaluationCallback : public ceres::EvaluationCallback {
+ public:
+  // m and c are passed by reference so that we have access to their values as
+  // they evolve over time through the course of optimization.
+  MyEvaluationCallback(const double& m, const double& c) : m_(m), c_(c) {
+    x_ = Eigen::VectorXd::Zero(kNumObservations);
+    y_ = Eigen::VectorXd::Zero(kNumObservations);
+    residuals_ = Eigen::VectorXd::Zero(kNumObservations);
+    jacobians_ = Eigen::MatrixXd::Zero(kNumObservations, 2);
+    for (int i = 0; i < kNumObservations; ++i) {
+      x_[i] = data[2 * i];
+      y_[i] = data[2 * i + 1];
+    }
+    PrepareForEvaluation(true, true);
+  }
+
+  void PrepareForEvaluation(bool evaluate_jacobians,
+                            bool new_evaluation_point) final {
+    if (new_evaluation_point) {
+      ComputeResidualAndJacobian(evaluate_jacobians);
+      jacobians_are_stale_ = !evaluate_jacobians;
+    } else {
+      if (evaluate_jacobians && jacobians_are_stale_) {
+        ComputeResidualAndJacobian(evaluate_jacobians);
+        jacobians_are_stale_ = false;
+      }
+    }
+  }
+
+  const Eigen::VectorXd& residuals() const { return residuals_; }
+  const Eigen::MatrixXd& jacobians() const { return jacobians_; }
+  bool jacobians_are_stale() const { return jacobians_are_stale_; }
+
+ private:
+  void ComputeResidualAndJacobian(bool evaluate_jacobians) {
+    residuals_ = -(m_ * x_.array() + c_).exp();
+    if (evaluate_jacobians) {
+      jacobians_.col(0) = residuals_.array() * x_.array();
+      jacobians_.col(1) = residuals_;
+    }
+    residuals_ += y_;
+  }
+
+  const double& m_;
+  const double& c_;
+  Eigen::VectorXd x_;
+  Eigen::VectorXd y_;
+  Eigen::VectorXd residuals_;
+  Eigen::MatrixXd jacobians_;
+
+  // jacobians_are_stale_ keeps track of whether the jacobian matrix matches the
+  // residuals or not, we only compute it if we know that Solver is going to
+  // need access to it.
+  bool jacobians_are_stale_ = true;
+};
+
+// As the name implies this CostFunction does not do any computation, it just
+// copies the appropriate residual and Jacobian from the matrices stored in
+// MyEvaluationCallback.
+class CostAndJacobianCopyingCostFunction
+    : public ceres::SizedCostFunction<1, 1, 1> {
+ public:
+  CostAndJacobianCopyingCostFunction(
+      int index, const MyEvaluationCallback& evaluation_callback)
+      : index_(index), evaluation_callback_(evaluation_callback) {}
+  ~CostAndJacobianCopyingCostFunction() override = default;
+
+  bool Evaluate(double const* const* parameters,
+                double* residuals,
+                double** jacobians) const final {
+    residuals[0] = evaluation_callback_.residuals()(index_);
+    if (!jacobians) return true;
+
+    // Ensure that we are not using stale Jacobians.
+    CHECK(!evaluation_callback_.jacobians_are_stale());
+
+    if (jacobians[0] != nullptr)
+      jacobians[0][0] = evaluation_callback_.jacobians()(index_, 0);
+    if (jacobians[1] != nullptr)
+      jacobians[1][0] = evaluation_callback_.jacobians()(index_, 1);
+    return true;
+  }
+
+ private:
+  int index_ = -1;
+  const MyEvaluationCallback& evaluation_callback_;
+};
+
+int main(int argc, char** argv) {
+  google::InitGoogleLogging(argv[0]);
+
+  const double initial_m = 0.0;
+  const double initial_c = 0.0;
+  double m = initial_m;
+  double c = initial_c;
+
+  MyEvaluationCallback evaluation_callback(m, c);
+  ceres::Problem::Options problem_options;
+  problem_options.evaluation_callback = &evaluation_callback;
+  ceres::Problem problem(problem_options);
+  for (int i = 0; i < kNumObservations; ++i) {
+    problem.AddResidualBlock(
+        new CostAndJacobianCopyingCostFunction(i, evaluation_callback),
+        nullptr,
+        &m,
+        &c);
+  }
+
+  ceres::Solver::Options options;
+  options.max_num_iterations = 25;
+  options.linear_solver_type = ceres::DENSE_QR;
+  options.minimizer_progress_to_stdout = true;
+
+  ceres::Solver::Summary summary;
+  ceres::Solve(options, &problem, &summary);
+  std::cout << summary.BriefReport() << "\n";
+  std::cout << "Initial m: " << initial_m << " c: " << initial_c << "\n";
+  std::cout << "Final   m: " << m << " c: " << c << "\n";
+  return 0;
+}
diff --git a/examples/fields_of_experts.cc b/examples/fields_of_experts.cc
index 7b7983e..f59fe16 100644
--- a/examples/fields_of_experts.cc
+++ b/examples/fields_of_experts.cc
@@ -1,5 +1,5 @@
 // Ceres Solver - A fast non-linear least squares minimizer
-// Copyright 2015 Google Inc. All rights reserved.
+// Copyright 2023 Google Inc. All rights reserved.
 // http://ceres-solver.org/
 //
 // Redistribution and use in source and binary forms, with or without
@@ -28,7 +28,7 @@
 //
 // Author: strandmark@google.com (Petter Strandmark)
 //
-// Class for loading the data required for descibing a Fields of Experts (FoE)
+// Class for loading the data required for describing a Fields of Experts (FoE)
 // model.
 
 #include "fields_of_experts.h"
@@ -38,13 +38,12 @@
 
 #include "pgm_image.h"
 
-namespace ceres {
-namespace examples {
+namespace ceres::examples {
 
 FieldsOfExpertsCost::FieldsOfExpertsCost(const std::vector<double>& filter)
     : filter_(filter) {
   set_num_residuals(1);
-  for (int i = 0; i < filter_.size(); ++i) {
+  for (int64_t i = 0; i < filter_.size(); ++i) {
     mutable_parameter_block_sizes()->push_back(1);
   }
 }
@@ -54,15 +53,15 @@
 bool FieldsOfExpertsCost::Evaluate(double const* const* parameters,
                                    double* residuals,
                                    double** jacobians) const {
-  int num_variables = filter_.size();
+  const int64_t num_variables = filter_.size();
   residuals[0] = 0;
-  for (int i = 0; i < num_variables; ++i) {
+  for (int64_t i = 0; i < num_variables; ++i) {
     residuals[0] += filter_[i] * parameters[i][0];
   }
 
-  if (jacobians != NULL) {
-    for (int i = 0; i < num_variables; ++i) {
-      if (jacobians[i] != NULL) {
+  if (jacobians != nullptr) {
+    for (int64_t i = 0; i < num_variables; ++i) {
+      if (jacobians[i] != nullptr) {
         jacobians[i][0] = filter_[i];
       }
     }
@@ -145,5 +144,4 @@
   return new FieldsOfExpertsLoss(alpha_[alpha_index]);
 }
 
-}  // namespace examples
-}  // namespace ceres
+}  // namespace ceres::examples
diff --git a/examples/fields_of_experts.h b/examples/fields_of_experts.h
index 429881d..2ff8c94 100644
--- a/examples/fields_of_experts.h
+++ b/examples/fields_of_experts.h
@@ -1,5 +1,5 @@
 // Ceres Solver - A fast non-linear least squares minimizer
-// Copyright 2015 Google Inc. All rights reserved.
+// Copyright 2023 Google Inc. All rights reserved.
 // http://ceres-solver.org/
 //
 // Redistribution and use in source and binary forms, with or without
@@ -28,7 +28,7 @@
 //
 // Author: strandmark@google.com (Petter Strandmark)
 //
-// Class for loading the data required for descibing a Fields of Experts (FoE)
+// Class for loading the data required for describing a Fields of Experts (FoE)
 // model. The Fields of Experts regularization consists of terms of the type
 //
 //   alpha * log(1 + (1/2)*sum(F .* X)^2),
@@ -52,8 +52,7 @@
 #include "ceres/sized_cost_function.h"
 #include "pgm_image.h"
 
-namespace ceres {
-namespace examples {
+namespace ceres::examples {
 
 // One sum in the FoE regularizer. This is a dot product between a filter and an
 // image patch. It simply calculates the dot product between the filter
@@ -63,9 +62,9 @@
   explicit FieldsOfExpertsCost(const std::vector<double>& filter);
   // The number of scalar parameters passed to Evaluate must equal the number of
   // filter coefficients passed to the constructor.
-  virtual bool Evaluate(double const* const* parameters,
-                        double* residuals,
-                        double** jacobians) const;
+  bool Evaluate(double const* const* parameters,
+                double* residuals,
+                double** jacobians) const override;
 
  private:
   const std::vector<double>& filter_;
@@ -78,7 +77,7 @@
 class FieldsOfExpertsLoss : public ceres::LossFunction {
  public:
   explicit FieldsOfExpertsLoss(double alpha) : alpha_(alpha) {}
-  virtual void Evaluate(double, double*) const;
+  void Evaluate(double, double*) const override;
 
  private:
   const double alpha_;
@@ -128,7 +127,6 @@
   std::vector<std::vector<double>> filters_;
 };
 
-}  // namespace examples
-}  // namespace ceres
+}  // namespace ceres::examples
 
 #endif  // CERES_EXAMPLES_FIELDS_OF_EXPERTS_H_
diff --git a/examples/helloworld.cc b/examples/helloworld.cc
index 9e32fad..40c2f2c 100644
--- a/examples/helloworld.cc
+++ b/examples/helloworld.cc
@@ -1,5 +1,5 @@
 // Ceres Solver - A fast non-linear least squares minimizer
-// Copyright 2015 Google Inc. All rights reserved.
+// Copyright 2023 Google Inc. All rights reserved.
 // http://ceres-solver.org/
 //
 // Redistribution and use in source and binary forms, with or without
@@ -36,12 +36,6 @@
 #include "ceres/ceres.h"
 #include "glog/logging.h"
 
-using ceres::AutoDiffCostFunction;
-using ceres::CostFunction;
-using ceres::Problem;
-using ceres::Solve;
-using ceres::Solver;
-
 // A templated cost functor that implements the residual r = 10 -
 // x. The method operator() is templated so that we can then use an
 // automatic differentiation wrapper around it to generate its
@@ -63,19 +57,19 @@
   const double initial_x = x;
 
   // Build the problem.
-  Problem problem;
+  ceres::Problem problem;
 
   // Set up the only cost function (also known as residual). This uses
   // auto-differentiation to obtain the derivative (jacobian).
-  CostFunction* cost_function =
-      new AutoDiffCostFunction<CostFunctor, 1, 1>(new CostFunctor);
+  ceres::CostFunction* cost_function =
+      new ceres::AutoDiffCostFunction<CostFunctor, 1, 1>();
   problem.AddResidualBlock(cost_function, nullptr, &x);
 
   // Run the solver!
-  Solver::Options options;
+  ceres::Solver::Options options;
   options.minimizer_progress_to_stdout = true;
-  Solver::Summary summary;
-  Solve(options, &problem, &summary);
+  ceres::Solver::Summary summary;
+  ceres::Solve(options, &problem, &summary);
 
   std::cout << summary.BriefReport() << "\n";
   std::cout << "x : " << initial_x << " -> " << x << "\n";
diff --git a/examples/helloworld_analytic_diff.cc b/examples/helloworld_analytic_diff.cc
index 6e120b5..b4826a2 100644
--- a/examples/helloworld_analytic_diff.cc
+++ b/examples/helloworld_analytic_diff.cc
@@ -1,5 +1,5 @@
 // Ceres Solver - A fast non-linear least squares minimizer
-// Copyright 2015 Google Inc. All rights reserved.
+// Copyright 2023 Google Inc. All rights reserved.
 // http://ceres-solver.org/
 //
 // Redistribution and use in source and binary forms, with or without
@@ -37,23 +37,15 @@
 #include "ceres/ceres.h"
 #include "glog/logging.h"
 
-using ceres::CostFunction;
-using ceres::Problem;
-using ceres::SizedCostFunction;
-using ceres::Solve;
-using ceres::Solver;
-
 // A CostFunction implementing analytically derivatives for the
 // function f(x) = 10 - x.
 class QuadraticCostFunction
-    : public SizedCostFunction<1 /* number of residuals */,
-                               1 /* size of first parameter */> {
+    : public ceres::SizedCostFunction<1 /* number of residuals */,
+                                      1 /* size of first parameter */> {
  public:
-  virtual ~QuadraticCostFunction() {}
-
-  virtual bool Evaluate(double const* const* parameters,
-                        double* residuals,
-                        double** jacobians) const {
+  bool Evaluate(double const* const* parameters,
+                double* residuals,
+                double** jacobians) const override {
     double x = parameters[0][0];
 
     // f(x) = 10 - x.
@@ -64,14 +56,14 @@
     // jacobians.
     //
     // Since the Evaluate function can be called with the jacobians
-    // pointer equal to NULL, the Evaluate function must check to see
+    // pointer equal to nullptr, the Evaluate function must check to see
     // if jacobians need to be computed.
     //
     // For this simple problem it is overkill to check if jacobians[0]
-    // is NULL, but in general when writing more complex
+    // is nullptr, but in general when writing more complex
     // CostFunctions, it is possible that Ceres may only demand the
     // derivatives w.r.t. a subset of the parameter blocks.
-    if (jacobians != NULL && jacobians[0] != NULL) {
+    if (jacobians != nullptr && jacobians[0] != nullptr) {
       jacobians[0][0] = -1;
     }
 
@@ -88,17 +80,17 @@
   const double initial_x = x;
 
   // Build the problem.
-  Problem problem;
+  ceres::Problem problem;
 
   // Set up the only cost function (also known as residual).
-  CostFunction* cost_function = new QuadraticCostFunction;
-  problem.AddResidualBlock(cost_function, NULL, &x);
+  ceres::CostFunction* cost_function = new QuadraticCostFunction;
+  problem.AddResidualBlock(cost_function, nullptr, &x);
 
   // Run the solver!
-  Solver::Options options;
+  ceres::Solver::Options options;
   options.minimizer_progress_to_stdout = true;
-  Solver::Summary summary;
-  Solve(options, &problem, &summary);
+  ceres::Solver::Summary summary;
+  ceres::Solve(options, &problem, &summary);
 
   std::cout << summary.BriefReport() << "\n";
   std::cout << "x : " << initial_x << " -> " << x << "\n";
diff --git a/examples/helloworld_numeric_diff.cc b/examples/helloworld_numeric_diff.cc
index 474adf3..4ed9ca6 100644
--- a/examples/helloworld_numeric_diff.cc
+++ b/examples/helloworld_numeric_diff.cc
@@ -1,5 +1,5 @@
 // Ceres Solver - A fast non-linear least squares minimizer
-// Copyright 2015 Google Inc. All rights reserved.
+// Copyright 2023 Google Inc. All rights reserved.
 // http://ceres-solver.org/
 //
 // Redistribution and use in source and binary forms, with or without
@@ -34,13 +34,6 @@
 #include "ceres/ceres.h"
 #include "glog/logging.h"
 
-using ceres::CENTRAL;
-using ceres::CostFunction;
-using ceres::NumericDiffCostFunction;
-using ceres::Problem;
-using ceres::Solve;
-using ceres::Solver;
-
 // A cost functor that implements the residual r = 10 - x.
 struct CostFunctor {
   bool operator()(const double* const x, double* residual) const {
@@ -58,19 +51,20 @@
   const double initial_x = x;
 
   // Build the problem.
-  Problem problem;
+  ceres::Problem problem;
 
   // Set up the only cost function (also known as residual). This uses
   // numeric differentiation to obtain the derivative (jacobian).
-  CostFunction* cost_function =
-      new NumericDiffCostFunction<CostFunctor, CENTRAL, 1, 1>(new CostFunctor);
-  problem.AddResidualBlock(cost_function, NULL, &x);
+  ceres::CostFunction* cost_function =
+      new ceres::NumericDiffCostFunction<CostFunctor, ceres::CENTRAL, 1, 1>(
+          new CostFunctor);
+  problem.AddResidualBlock(cost_function, nullptr, &x);
 
   // Run the solver!
-  Solver::Options options;
+  ceres::Solver::Options options;
   options.minimizer_progress_to_stdout = true;
-  Solver::Summary summary;
-  Solve(options, &problem, &summary);
+  ceres::Solver::Summary summary;
+  ceres::Solve(options, &problem, &summary);
 
   std::cout << summary.BriefReport() << "\n";
   std::cout << "x : " << initial_x << " -> " << x << "\n";
diff --git a/examples/iteration_callback_example.cc b/examples/iteration_callback_example.cc
new file mode 100644
index 0000000..0be2f36
--- /dev/null
+++ b/examples/iteration_callback_example.cc
@@ -0,0 +1,199 @@
+// Ceres Solver - A fast non-linear least squares minimizer
+// Copyright 2023 Google Inc. All rights reserved.
+// http://ceres-solver.org/
+//
+// Redistribution and use in source and binary forms, with or without
+// modification, are permitted provided that the following conditions are met:
+//
+// * Redistributions of source code must retain the above copyright notice,
+//   this list of conditions and the following disclaimer.
+// * Redistributions in binary form must reproduce the above copyright notice,
+//   this list of conditions and the following disclaimer in the documentation
+//   and/or other materials provided with the distribution.
+// * Neither the name of Google Inc. nor the names of its contributors may be
+//   used to endorse or promote products derived from this software without
+//   specific prior written permission.
+//
+// THIS SOFTWARE IS PROVIDED BY THE COPYRIGHT HOLDERS AND CONTRIBUTORS "AS IS"
+// AND ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED TO, THE
+// IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR PURPOSE
+// ARE DISCLAIMED. IN NO EVENT SHALL THE COPYRIGHT OWNER OR CONTRIBUTORS BE
+// LIABLE FOR ANY DIRECT, INDIRECT, INCIDENTAL, SPECIAL, EXEMPLARY, OR
+// CONSEQUENTIAL DAMAGES (INCLUDING, BUT NOT LIMITED TO, PROCUREMENT OF
+// SUBSTITUTE GOODS OR SERVICES; LOSS OF USE, DATA, OR PROFITS; OR BUSINESS
+// INTERRUPTION) HOWEVER CAUSED AND ON ANY THEORY OF LIABILITY, WHETHER IN
+// CONTRACT, STRICT LIABILITY, OR TORT (INCLUDING NEGLIGENCE OR OTHERWISE)
+// ARISING IN ANY WAY OUT OF THE USE OF THIS SOFTWARE, EVEN IF ADVISED OF THE
+// POSSIBILITY OF SUCH DAMAGE.
+//
+// Author: sameeragarwal@google.com (Sameer Agarwal)
+//
+// This example is a variant of curve_fitting.cc where we use an
+// IterationCallback to implement custom logging which prints out the values of
+// the parameter blocks as they evolve over the course of the optimization. This
+// also requires the use of Solver::Options::update_state_every_iteration.
+
+#include <iostream>
+
+#include "ceres/ceres.h"
+#include "glog/logging.h"
+
+// Data generated using the following octave code.
+//   randn('seed', 23497);
+//   m = 0.3;
+//   c = 0.1;
+//   x=[0:0.075:5];
+//   y = exp(m * x + c);
+//   noise = randn(size(x)) * 0.2;
+//   y_observed = y + noise;
+//   data = [x', y_observed'];
+
+const int kNumObservations = 67;
+// clang-format off
+const double data[] = {
+  0.000000e+00, 1.133898e+00,
+  7.500000e-02, 1.334902e+00,
+  1.500000e-01, 1.213546e+00,
+  2.250000e-01, 1.252016e+00,
+  3.000000e-01, 1.392265e+00,
+  3.750000e-01, 1.314458e+00,
+  4.500000e-01, 1.472541e+00,
+  5.250000e-01, 1.536218e+00,
+  6.000000e-01, 1.355679e+00,
+  6.750000e-01, 1.463566e+00,
+  7.500000e-01, 1.490201e+00,
+  8.250000e-01, 1.658699e+00,
+  9.000000e-01, 1.067574e+00,
+  9.750000e-01, 1.464629e+00,
+  1.050000e+00, 1.402653e+00,
+  1.125000e+00, 1.713141e+00,
+  1.200000e+00, 1.527021e+00,
+  1.275000e+00, 1.702632e+00,
+  1.350000e+00, 1.423899e+00,
+  1.425000e+00, 1.543078e+00,
+  1.500000e+00, 1.664015e+00,
+  1.575000e+00, 1.732484e+00,
+  1.650000e+00, 1.543296e+00,
+  1.725000e+00, 1.959523e+00,
+  1.800000e+00, 1.685132e+00,
+  1.875000e+00, 1.951791e+00,
+  1.950000e+00, 2.095346e+00,
+  2.025000e+00, 2.361460e+00,
+  2.100000e+00, 2.169119e+00,
+  2.175000e+00, 2.061745e+00,
+  2.250000e+00, 2.178641e+00,
+  2.325000e+00, 2.104346e+00,
+  2.400000e+00, 2.584470e+00,
+  2.475000e+00, 1.914158e+00,
+  2.550000e+00, 2.368375e+00,
+  2.625000e+00, 2.686125e+00,
+  2.700000e+00, 2.712395e+00,
+  2.775000e+00, 2.499511e+00,
+  2.850000e+00, 2.558897e+00,
+  2.925000e+00, 2.309154e+00,
+  3.000000e+00, 2.869503e+00,
+  3.075000e+00, 3.116645e+00,
+  3.150000e+00, 3.094907e+00,
+  3.225000e+00, 2.471759e+00,
+  3.300000e+00, 3.017131e+00,
+  3.375000e+00, 3.232381e+00,
+  3.450000e+00, 2.944596e+00,
+  3.525000e+00, 3.385343e+00,
+  3.600000e+00, 3.199826e+00,
+  3.675000e+00, 3.423039e+00,
+  3.750000e+00, 3.621552e+00,
+  3.825000e+00, 3.559255e+00,
+  3.900000e+00, 3.530713e+00,
+  3.975000e+00, 3.561766e+00,
+  4.050000e+00, 3.544574e+00,
+  4.125000e+00, 3.867945e+00,
+  4.200000e+00, 4.049776e+00,
+  4.275000e+00, 3.885601e+00,
+  4.350000e+00, 4.110505e+00,
+  4.425000e+00, 4.345320e+00,
+  4.500000e+00, 4.161241e+00,
+  4.575000e+00, 4.363407e+00,
+  4.650000e+00, 4.161576e+00,
+  4.725000e+00, 4.619728e+00,
+  4.800000e+00, 4.737410e+00,
+  4.875000e+00, 4.727863e+00,
+  4.950000e+00, 4.669206e+00,
+};
+// clang-format on
+
+struct ExponentialResidual {
+  ExponentialResidual(double x, double y) : x(x), y(y) {}
+
+  template <typename T>
+  bool operator()(const T* const m, const T* const c, T* residual) const {
+    residual[0] = y - exp(m[0] * x + c[0]);
+    return true;
+  }
+
+ private:
+  const double x;
+  const double y;
+};
+
+// MyIterationCallback prints the iteration number, the cost and the value of
+// the parameter blocks every iteration.
+class MyIterationCallback : public ceres::IterationCallback {
+ public:
+  MyIterationCallback(const double* m, const double* c) : m_(m), c_(c) {}
+
+  ~MyIterationCallback() override = default;
+
+  ceres::CallbackReturnType operator()(
+      const ceres::IterationSummary& summary) final {
+    std::cout << "Iteration: " << summary.iteration << " cost: " << summary.cost
+              << " m: " << *m_ << " c: " << *c_ << std::endl;
+    return ceres::SOLVER_CONTINUE;
+  }
+
+ private:
+  const double* m_ = nullptr;
+  const double* c_ = nullptr;
+};
+
+int main(int argc, char** argv) {
+  google::InitGoogleLogging(argv[0]);
+
+  const double initial_m = 0.0;
+  const double initial_c = 0.0;
+
+  double m = initial_m;
+  double c = initial_c;
+
+  ceres::Problem problem;
+  for (int i = 0; i < kNumObservations; ++i) {
+    problem.AddResidualBlock(
+        new ceres::AutoDiffCostFunction<ExponentialResidual, 1, 1, 1>(
+            data[2 * i], data[2 * i + 1]),
+        nullptr,
+        &m,
+        &c);
+  }
+
+  ceres::Solver::Options options;
+  options.max_num_iterations = 25;
+  options.linear_solver_type = ceres::DENSE_QR;
+
+  // Turn off the default logging from Ceres so that it does not interfere with
+  // MyIterationCallback.
+  options.minimizer_progress_to_stdout = false;
+
+  MyIterationCallback callback(&m, &c);
+  options.callbacks.push_back(&callback);
+
+  // Tell Ceres to update the value of the parameter blocks on each each
+  // iteration (successful or not) so that MyIterationCallback will be able to
+  // see them when called.
+  options.update_state_every_iteration = true;
+
+  ceres::Solver::Summary summary;
+  ceres::Solve(options, &problem, &summary);
+  std::cout << summary.BriefReport() << "\n";
+  std::cout << "Initial m: " << initial_m << " c: " << initial_c << "\n";
+  std::cout << "Final   m: " << m << " c: " << c << "\n";
+  return 0;
+}
diff --git a/examples/libmv_bundle_adjuster.cc b/examples/libmv_bundle_adjuster.cc
index b1eb220..9315ed7 100644
--- a/examples/libmv_bundle_adjuster.cc
+++ b/examples/libmv_bundle_adjuster.cc
@@ -60,7 +60,7 @@
 // Image number shall be greater or equal to zero. Order of cameras does not
 // matter and gaps are possible.
 //
-// Every 3D point is decribed by:
+// Every 3D point is described by:
 //
 //  - Track number point belongs to (single 4 bytes integer value).
 //  - 3D position vector, 3-component vector of float values.
@@ -100,12 +100,16 @@
 #define close _close
 typedef unsigned __int32 uint32_t;
 #else
-#include <stdint.h>
 #include <unistd.h>
 
+#include <cstdint>
+
+// NOTE MinGW does define the macro.
+#ifndef O_BINARY
 // O_BINARY is not defined on unix like platforms, as there is no
 // difference between binary and text files.
 #define O_BINARY 0
+#endif
 
 #endif
 
@@ -114,12 +118,10 @@
 #include "gflags/gflags.h"
 #include "glog/logging.h"
 
-typedef Eigen::Matrix<double, 3, 3> Mat3;
-typedef Eigen::Matrix<double, 6, 1> Vec6;
-typedef Eigen::Vector3d Vec3;
-typedef Eigen::Vector4d Vec4;
-
-using std::vector;
+using Mat3 = Eigen::Matrix<double, 3, 3>;
+using Vec6 = Eigen::Matrix<double, 6, 1>;
+using Vec3 = Eigen::Vector3d;
+using Vec4 = Eigen::Vector4d;
 
 DEFINE_string(input, "", "Input File name");
 DEFINE_string(refine_intrinsics,
@@ -135,10 +137,10 @@
 // R is a 3x3 matrix representing the rotation of the camera.
 // t is a translation vector representing its positions.
 struct EuclideanCamera {
-  EuclideanCamera() : image(-1) {}
-  EuclideanCamera(const EuclideanCamera& c) : image(c.image), R(c.R), t(c.t) {}
+  EuclideanCamera() = default;
+  EuclideanCamera(const EuclideanCamera& c) = default;
 
-  int image;
+  int image{-1};
   Mat3 R;
   Vec3 t;
 };
@@ -148,9 +150,9 @@
 // track identifies which track this point corresponds to.
 // X represents the 3D position of the track.
 struct EuclideanPoint {
-  EuclideanPoint() : track(-1) {}
-  EuclideanPoint(const EuclideanPoint& p) : track(p.track), X(p.X) {}
-  int track;
+  EuclideanPoint() = default;
+  EuclideanPoint(const EuclideanPoint& p) = default;
+  int track{-1};
   Vec3 X;
 };
 
@@ -203,32 +205,32 @@
 };
 
 // Returns a pointer to the camera corresponding to a image.
-EuclideanCamera* CameraForImage(vector<EuclideanCamera>* all_cameras,
+EuclideanCamera* CameraForImage(std::vector<EuclideanCamera>* all_cameras,
                                 const int image) {
   if (image < 0 || image >= all_cameras->size()) {
-    return NULL;
+    return nullptr;
   }
   EuclideanCamera* camera = &(*all_cameras)[image];
   if (camera->image == -1) {
-    return NULL;
+    return nullptr;
   }
   return camera;
 }
 
 const EuclideanCamera* CameraForImage(
-    const vector<EuclideanCamera>& all_cameras, const int image) {
+    const std::vector<EuclideanCamera>& all_cameras, const int image) {
   if (image < 0 || image >= all_cameras.size()) {
-    return NULL;
+    return nullptr;
   }
   const EuclideanCamera* camera = &all_cameras[image];
   if (camera->image == -1) {
-    return NULL;
+    return nullptr;
   }
   return camera;
 }
 
 // Returns maximal image number at which marker exists.
-int MaxImage(const vector<Marker>& all_markers) {
+int MaxImage(const std::vector<Marker>& all_markers) {
   if (all_markers.size() == 0) {
     return -1;
   }
@@ -241,14 +243,14 @@
 }
 
 // Returns a pointer to the point corresponding to a track.
-EuclideanPoint* PointForTrack(vector<EuclideanPoint>* all_points,
+EuclideanPoint* PointForTrack(std::vector<EuclideanPoint>* all_points,
                               const int track) {
   if (track < 0 || track >= all_points->size()) {
-    return NULL;
+    return nullptr;
   }
   EuclideanPoint* point = &(*all_points)[track];
   if (point->track == -1) {
-    return NULL;
+    return nullptr;
   }
   return point;
 }
@@ -262,7 +264,7 @@
 // denotes file endianness in this way.
 class EndianAwareFileReader {
  public:
-  EndianAwareFileReader(void) : file_descriptor_(-1) {
+  EndianAwareFileReader() {
     // Get an endian type of the host machine.
     union {
       unsigned char bytes[4];
@@ -272,7 +274,7 @@
     file_endian_type_ = host_endian_type_;
   }
 
-  ~EndianAwareFileReader(void) {
+  ~EndianAwareFileReader() {
     if (file_descriptor_ > 0) {
       close(file_descriptor_);
     }
@@ -284,7 +286,7 @@
       return false;
     }
     // Get an endian tpye of data in the file.
-    unsigned char file_endian_type_flag = Read<unsigned char>();
+    auto file_endian_type_flag = Read<unsigned char>();
     if (file_endian_type_flag == 'V') {
       file_endian_type_ = kBigEndian;
     } else if (file_endian_type_flag == 'v') {
@@ -297,9 +299,11 @@
 
   // Read value from the file, will switch endian if needed.
   template <typename T>
-  T Read(void) const {
+  T Read() const {
     T value;
+    CERES_DISABLE_DEPRECATED_WARNING
     CHECK_GT(read(file_descriptor_, &value, sizeof(value)), 0);
+    CERES_RESTORE_DEPRECATED_WARNING
     // Switch endian type if file contains data in different type
     // that current machine.
     if (file_endian_type_ != host_endian_type_) {
@@ -316,7 +320,7 @@
   template <typename T>
   T SwitchEndian(const T value) const {
     if (sizeof(T) == 4) {
-      unsigned int temp_value = static_cast<unsigned int>(value);
+      auto temp_value = static_cast<unsigned int>(value);
       // clang-format off
       return ((temp_value >> 24)) |
              ((temp_value << 8) & 0x00ff0000) |
@@ -333,7 +337,7 @@
 
   int host_endian_type_;
   int file_endian_type_;
-  int file_descriptor_;
+  int file_descriptor_{-1};
 };
 
 // Read 3x3 column-major matrix from the file
@@ -369,17 +373,17 @@
 // reading.
 bool ReadProblemFromFile(const std::string& file_name,
                          double camera_intrinsics[8],
-                         vector<EuclideanCamera>* all_cameras,
-                         vector<EuclideanPoint>* all_points,
+                         std::vector<EuclideanCamera>* all_cameras,
+                         std::vector<EuclideanPoint>* all_points,
                          bool* is_image_space,
-                         vector<Marker>* all_markers) {
+                         std::vector<Marker>* all_markers) {
   EndianAwareFileReader file_reader;
   if (!file_reader.OpenFile(file_name)) {
     return false;
   }
 
   // Read markers' space flag.
-  unsigned char is_image_space_flag = file_reader.Read<unsigned char>();
+  auto is_image_space_flag = file_reader.Read<unsigned char>();
   if (is_image_space_flag == 'P') {
     *is_image_space = true;
   } else if (is_image_space_flag == 'N') {
@@ -610,10 +614,10 @@
 //
 // Element with index i matches to a rotation+translation for
 // camera at image i.
-vector<Vec6> PackCamerasRotationAndTranslation(
-    const vector<Marker>& all_markers,
-    const vector<EuclideanCamera>& all_cameras) {
-  vector<Vec6> all_cameras_R_t;
+std::vector<Vec6> PackCamerasRotationAndTranslation(
+    const std::vector<Marker>& all_markers,
+    const std::vector<EuclideanCamera>& all_cameras) {
+  std::vector<Vec6> all_cameras_R_t;
   int max_image = MaxImage(all_markers);
 
   all_cameras_R_t.resize(max_image + 1);
@@ -633,9 +637,10 @@
 }
 
 // Convert cameras rotations fro mangle axis back to rotation matrix.
-void UnpackCamerasRotationAndTranslation(const vector<Marker>& all_markers,
-                                         const vector<Vec6>& all_cameras_R_t,
-                                         vector<EuclideanCamera>* all_cameras) {
+void UnpackCamerasRotationAndTranslation(
+    const std::vector<Marker>& all_markers,
+    const std::vector<Vec6>& all_cameras_R_t,
+    std::vector<EuclideanCamera>* all_cameras) {
   int max_image = MaxImage(all_markers);
 
   for (int i = 0; i <= max_image; i++) {
@@ -650,12 +655,12 @@
   }
 }
 
-void EuclideanBundleCommonIntrinsics(const vector<Marker>& all_markers,
+void EuclideanBundleCommonIntrinsics(const std::vector<Marker>& all_markers,
                                      const int bundle_intrinsics,
                                      const int bundle_constraints,
                                      double* camera_intrinsics,
-                                     vector<EuclideanCamera>* all_cameras,
-                                     vector<EuclideanPoint>* all_points) {
+                                     std::vector<EuclideanCamera>* all_cameras,
+                                     std::vector<EuclideanPoint>* all_points) {
   PrintCameraIntrinsics("Original intrinsics: ", camera_intrinsics);
 
   ceres::Problem::Options problem_options;
@@ -667,11 +672,11 @@
   //
   // Block for minimization has got the following structure:
   //   <3 elements for angle-axis> <3 elements for translation>
-  vector<Vec6> all_cameras_R_t =
+  std::vector<Vec6> all_cameras_R_t =
       PackCamerasRotationAndTranslation(all_markers, *all_cameras);
 
-  // Parameterization used to restrict camera motion for modal solvers.
-  ceres::SubsetParameterization* constant_transform_parameterization = NULL;
+  // Manifold used to restrict camera motion for modal solvers.
+  ceres::SubsetManifold* constant_transform_manifold = nullptr;
   if (bundle_constraints & BUNDLE_NO_TRANSLATION) {
     std::vector<int> constant_translation;
 
@@ -680,8 +685,8 @@
     constant_translation.push_back(4);
     constant_translation.push_back(5);
 
-    constant_transform_parameterization =
-        new ceres::SubsetParameterization(6, constant_translation);
+    constant_transform_manifold =
+        new ceres::SubsetManifold(6, constant_translation);
   }
 
   std::vector<OpenCVReprojectionError> errors;
@@ -692,11 +697,10 @@
 
   int num_residuals = 0;
   bool have_locked_camera = false;
-  for (int i = 0; i < all_markers.size(); ++i) {
-    const Marker& marker = all_markers[i];
+  for (const auto& marker : all_markers) {
     EuclideanCamera* camera = CameraForImage(all_cameras, marker.image);
     EuclideanPoint* point = PointForTrack(all_points, marker.track);
-    if (camera == NULL || point == NULL) {
+    if (camera == nullptr || point == nullptr) {
       continue;
     }
 
@@ -708,7 +712,7 @@
     costFunctions.emplace_back(&errors.back(), ceres::DO_NOT_TAKE_OWNERSHIP);
 
     problem.AddResidualBlock(&costFunctions.back(),
-                             NULL,
+                             nullptr,
                              camera_intrinsics,
                              current_camera_R_t,
                              &point->X(0));
@@ -720,8 +724,7 @@
     }
 
     if (bundle_constraints & BUNDLE_NO_TRANSLATION) {
-      problem.SetParameterization(current_camera_R_t,
-                                  constant_transform_parameterization);
+      problem.SetManifold(current_camera_R_t, constant_transform_manifold);
     }
 
     num_residuals++;
@@ -760,10 +763,8 @@
     // Always set K3 constant, it's not used at the moment.
     constant_intrinsics.push_back(OFFSET_K3);
 
-    ceres::SubsetParameterization* subset_parameterization =
-        new ceres::SubsetParameterization(8, constant_intrinsics);
-
-    problem.SetParameterization(camera_intrinsics, subset_parameterization);
+    auto* subset_manifold = new ceres::SubsetManifold(8, constant_intrinsics);
+    problem.SetManifold(camera_intrinsics, subset_manifold);
   }
 
   // Configure the solver.
@@ -793,18 +794,18 @@
   GFLAGS_NAMESPACE::ParseCommandLineFlags(&argc, &argv, true);
   google::InitGoogleLogging(argv[0]);
 
-  if (FLAGS_input.empty()) {
+  if (CERES_GET_FLAG(FLAGS_input).empty()) {
     LOG(ERROR) << "Usage: libmv_bundle_adjuster --input=blender_problem";
     return EXIT_FAILURE;
   }
 
   double camera_intrinsics[8];
-  vector<EuclideanCamera> all_cameras;
-  vector<EuclideanPoint> all_points;
+  std::vector<EuclideanCamera> all_cameras;
+  std::vector<EuclideanPoint> all_points;
   bool is_image_space;
-  vector<Marker> all_markers;
+  std::vector<Marker> all_markers;
 
-  if (!ReadProblemFromFile(FLAGS_input,
+  if (!ReadProblemFromFile(CERES_GET_FLAG(FLAGS_input),
                            camera_intrinsics,
                            &all_cameras,
                            &all_points,
@@ -828,14 +829,14 @@
   // declare which intrinsics need to be refined and in this case
   // refining flags does not depend on problem at all.
   int bundle_intrinsics = BUNDLE_NO_INTRINSICS;
-  if (FLAGS_refine_intrinsics.empty()) {
+  if (CERES_GET_FLAG(FLAGS_refine_intrinsics).empty()) {
     if (is_image_space) {
       bundle_intrinsics = BUNDLE_FOCAL_LENGTH | BUNDLE_RADIAL;
     }
   } else {
-    if (FLAGS_refine_intrinsics == "radial") {
+    if (CERES_GET_FLAG(FLAGS_refine_intrinsics) == "radial") {
       bundle_intrinsics = BUNDLE_FOCAL_LENGTH | BUNDLE_RADIAL;
-    } else if (FLAGS_refine_intrinsics != "none") {
+    } else if (CERES_GET_FLAG(FLAGS_refine_intrinsics) != "none") {
       LOG(ERROR) << "Unsupported value for refine-intrinsics";
       return EXIT_FAILURE;
     }
diff --git a/examples/libmv_homography.cc b/examples/libmv_homography.cc
index 55f3b70..b7c9eda 100644
--- a/examples/libmv_homography.cc
+++ b/examples/libmv_homography.cc
@@ -1,5 +1,5 @@
 // Ceres Solver - A fast non-linear least squares minimizer
-// Copyright 2015 Google Inc. All rights reserved.
+// Copyright 2023 Google Inc. All rights reserved.
 // http://ceres-solver.org/
 //
 // Redistribution and use in source and binary forms, with or without
@@ -60,17 +60,19 @@
 // This example demonstrates custom exit criterion by having a callback check
 // for image-space error.
 
+#include <utility>
+
 #include "ceres/ceres.h"
 #include "glog/logging.h"
 
-typedef Eigen::NumTraits<double> EigenDouble;
+using EigenDouble = Eigen::NumTraits<double>;
 
-typedef Eigen::MatrixXd Mat;
-typedef Eigen::VectorXd Vec;
-typedef Eigen::Matrix<double, 3, 3> Mat3;
-typedef Eigen::Matrix<double, 2, 1> Vec2;
-typedef Eigen::Matrix<double, Eigen::Dynamic, 8> MatX8;
-typedef Eigen::Vector3d Vec3;
+using Mat = Eigen::MatrixXd;
+using Vec = Eigen::VectorXd;
+using Mat3 = Eigen::Matrix<double, 3, 3>;
+using Vec2 = Eigen::Matrix<double, 2, 1>;
+using MatX8 = Eigen::Matrix<double, Eigen::Dynamic, 8>;
+using Vec3 = Eigen::Vector3d;
 
 namespace {
 
@@ -82,11 +84,10 @@
 struct EstimateHomographyOptions {
   // Default settings for homography estimation which should be suitable
   // for a wide range of use cases.
-  EstimateHomographyOptions()
-      : max_num_iterations(50), expected_average_symmetric_distance(1e-16) {}
+  EstimateHomographyOptions() = default;
 
   // Maximal number of iterations for the refinement step.
-  int max_num_iterations;
+  int max_num_iterations{50};
 
   // Expected average of symmetric geometric distance between
   // actual destination points and original ones transformed by
@@ -96,7 +97,7 @@
   // geometric distance is less or equal to this value.
   //
   // This distance is measured in the same units as input points are.
-  double expected_average_symmetric_distance;
+  double expected_average_symmetric_distance{1e-16};
 };
 
 // Calculate symmetric geometric cost terms:
@@ -111,7 +112,7 @@
                                      const Eigen::Matrix<T, 2, 1>& x2,
                                      T forward_error[2],
                                      T backward_error[2]) {
-  typedef Eigen::Matrix<T, 3, 1> Vec3;
+  using Vec3 = Eigen::Matrix<T, 3, 1>;
   Vec3 x(x1(0), x1(1), T(1.0));
   Vec3 y(x2(0), x2(1), T(1.0));
 
@@ -152,8 +153,8 @@
 template <typename T = double>
 class Homography2DNormalizedParameterization {
  public:
-  typedef Eigen::Matrix<T, 8, 1> Parameters;     // a, b, ... g, h
-  typedef Eigen::Matrix<T, 3, 3> Parameterized;  // H
+  using Parameters = Eigen::Matrix<T, 8, 1>;     // a, b, ... g, h
+  using Parameterized = Eigen::Matrix<T, 3, 3>;  // H
 
   // Convert from the 8 parameters to a H matrix.
   static void To(const Parameters& p, Parameterized* h) {
@@ -202,11 +203,11 @@
   assert(x1.rows() == x2.rows());
   assert(x1.cols() == x2.cols());
 
-  int n = x1.cols();
+  const int64_t n = x1.cols();
   MatX8 L = Mat::Zero(n * 3, 8);
   Mat b = Mat::Zero(n * 3, 1);
-  for (int i = 0; i < n; ++i) {
-    int j = 3 * i;
+  for (int64_t i = 0; i < n; ++i) {
+    int64_t j = 3 * i;
     L(j, 0) = x1(0, i);              // a
     L(j, 1) = x1(1, i);              // b
     L(j, 2) = 1.0;                   // c
@@ -242,13 +243,13 @@
 // used for homography matrix refinement.
 class HomographySymmetricGeometricCostFunctor {
  public:
-  HomographySymmetricGeometricCostFunctor(const Vec2& x, const Vec2& y)
-      : x_(x), y_(y) {}
+  HomographySymmetricGeometricCostFunctor(Vec2 x, Vec2 y)
+      : x_(std::move(x)), y_(std::move(y)) {}
 
   template <typename T>
   bool operator()(const T* homography_parameters, T* residuals) const {
-    typedef Eigen::Matrix<T, 3, 3> Mat3;
-    typedef Eigen::Matrix<T, 2, 1> Vec2;
+    using Mat3 = Eigen::Matrix<T, 3, 3>;
+    using Vec2 = Eigen::Matrix<T, 2, 1>;
 
     Mat3 H(homography_parameters);
     Vec2 x(T(x_(0)), T(x_(1)));
@@ -277,8 +278,8 @@
                               Mat3* H)
       : options_(options), x1_(x1), x2_(x2), H_(H) {}
 
-  virtual ceres::CallbackReturnType operator()(
-      const ceres::IterationSummary& summary) {
+  ceres::CallbackReturnType operator()(
+      const ceres::IterationSummary& summary) override {
     // If the step wasn't successful, there's nothing to do.
     if (!summary.step_is_successful) {
       return ceres::SOLVER_CONTINUE;
@@ -326,16 +327,11 @@
   // Step 2: Refine matrix using Ceres minimizer.
   ceres::Problem problem;
   for (int i = 0; i < x1.cols(); i++) {
-    HomographySymmetricGeometricCostFunctor*
-        homography_symmetric_geometric_cost_function =
-            new HomographySymmetricGeometricCostFunctor(x1.col(i), x2.col(i));
-
     problem.AddResidualBlock(
         new ceres::AutoDiffCostFunction<HomographySymmetricGeometricCostFunctor,
                                         4,  // num_residuals
-                                        9>(
-            homography_symmetric_geometric_cost_function),
-        NULL,
+                                        9>(x1.col(i), x2.col(i)),
+        nullptr,
         H->data());
   }
 
@@ -380,10 +376,10 @@
 
   Mat x2 = x1;
   for (int i = 0; i < x2.cols(); ++i) {
-    Vec3 homogenous_x1 = Vec3(x1(0, i), x1(1, i), 1.0);
-    Vec3 homogenous_x2 = homography_matrix * homogenous_x1;
-    x2(0, i) = homogenous_x2(0) / homogenous_x2(2);
-    x2(1, i) = homogenous_x2(1) / homogenous_x2(2);
+    Vec3 homogeneous_x1 = Vec3(x1(0, i), x1(1, i), 1.0);
+    Vec3 homogeneous_x2 = homography_matrix * homogeneous_x1;
+    x2(0, i) = homogeneous_x2(0) / homogeneous_x2(2);
+    x2(1, i) = homogeneous_x2(1) / homogeneous_x2(2);
 
     // Apply some noise so algebraic estimation is not good enough.
     x2(0, i) += static_cast<double>(rand() % 1000) / 5000.0;
diff --git a/examples/more_garbow_hillstrom.cc b/examples/more_garbow_hillstrom.cc
index e39d23c..f15e576 100644
--- a/examples/more_garbow_hillstrom.cc
+++ b/examples/more_garbow_hillstrom.cc
@@ -1,5 +1,5 @@
 // Ceres Solver - A fast non-linear least squares minimizer
-// Copyright 2015 Google Inc. All rights reserved.
+// Copyright 2023 Google Inc. All rights reserved.
 // http://ceres-solver.org/
 //
 // Redistribution and use in source and binary forms, with or without
@@ -72,55 +72,56 @@
              3,
              "Maximal number of extrapolations in Ridders' method.");
 
-namespace ceres {
-namespace examples {
+namespace ceres::examples {
 
 const double kDoubleMax = std::numeric_limits<double>::max();
 
 static void SetNumericDiffOptions(ceres::NumericDiffOptions* options) {
-  options->max_num_ridders_extrapolations = FLAGS_ridders_extrapolations;
+  options->max_num_ridders_extrapolations =
+      CERES_GET_FLAG(FLAGS_ridders_extrapolations);
 }
 
-#define BEGIN_MGH_PROBLEM(name, num_parameters, num_residuals)                \
-  struct name {                                                               \
-    static constexpr int kNumParameters = num_parameters;                     \
-    static const double initial_x[kNumParameters];                            \
-    static const double lower_bounds[kNumParameters];                         \
-    static const double upper_bounds[kNumParameters];                         \
-    static const double constrained_optimal_cost;                             \
-    static const double unconstrained_optimal_cost;                           \
-    static CostFunction* Create() {                                           \
-      if (FLAGS_use_numeric_diff) {                                           \
-        ceres::NumericDiffOptions options;                                    \
-        SetNumericDiffOptions(&options);                                      \
-        if (FLAGS_numeric_diff_method == "central") {                         \
-          return new NumericDiffCostFunction<name,                            \
-                                             ceres::CENTRAL,                  \
-                                             num_residuals,                   \
-                                             num_parameters>(                 \
-              new name, ceres::TAKE_OWNERSHIP, num_residuals, options);       \
-        } else if (FLAGS_numeric_diff_method == "forward") {                  \
-          return new NumericDiffCostFunction<name,                            \
-                                             ceres::FORWARD,                  \
-                                             num_residuals,                   \
-                                             num_parameters>(                 \
-              new name, ceres::TAKE_OWNERSHIP, num_residuals, options);       \
-        } else if (FLAGS_numeric_diff_method == "ridders") {                  \
-          return new NumericDiffCostFunction<name,                            \
-                                             ceres::RIDDERS,                  \
-                                             num_residuals,                   \
-                                             num_parameters>(                 \
-              new name, ceres::TAKE_OWNERSHIP, num_residuals, options);       \
-        } else {                                                              \
-          LOG(ERROR) << "Invalid numeric diff method specified";              \
-          return NULL;                                                        \
-        }                                                                     \
-      } else {                                                                \
-        return new AutoDiffCostFunction<name, num_residuals, num_parameters>( \
-            new name);                                                        \
-      }                                                                       \
-    }                                                                         \
-    template <typename T>                                                     \
+#define BEGIN_MGH_PROBLEM(name, num_parameters, num_residuals)               \
+  struct name {                                                              \
+    static constexpr int kNumParameters = num_parameters;                    \
+    static const double initial_x[kNumParameters];                           \
+    static const double lower_bounds[kNumParameters];                        \
+    static const double upper_bounds[kNumParameters];                        \
+    static const double constrained_optimal_cost;                            \
+    static const double unconstrained_optimal_cost;                          \
+    static CostFunction* Create() {                                          \
+      if (CERES_GET_FLAG(FLAGS_use_numeric_diff)) {                          \
+        ceres::NumericDiffOptions options;                                   \
+        SetNumericDiffOptions(&options);                                     \
+        if (CERES_GET_FLAG(FLAGS_numeric_diff_method) == "central") {        \
+          return new NumericDiffCostFunction<name,                           \
+                                             ceres::CENTRAL,                 \
+                                             num_residuals,                  \
+                                             num_parameters>(                \
+              new name, ceres::TAKE_OWNERSHIP, num_residuals, options);      \
+        } else if (CERES_GET_FLAG(FLAGS_numeric_diff_method) == "forward") { \
+          return new NumericDiffCostFunction<name,                           \
+                                             ceres::FORWARD,                 \
+                                             num_residuals,                  \
+                                             num_parameters>(                \
+              new name, ceres::TAKE_OWNERSHIP, num_residuals, options);      \
+        } else if (CERES_GET_FLAG(FLAGS_numeric_diff_method) == "ridders") { \
+          return new NumericDiffCostFunction<name,                           \
+                                             ceres::RIDDERS,                 \
+                                             num_residuals,                  \
+                                             num_parameters>(                \
+              new name, ceres::TAKE_OWNERSHIP, num_residuals, options);      \
+        } else {                                                             \
+          LOG(ERROR) << "Invalid numeric diff method specified";             \
+          return nullptr;                                                    \
+        }                                                                    \
+      } else {                                                               \
+        return new AutoDiffCostFunction<name,                                \
+                                        num_residuals,                       \
+                                        num_parameters>();                   \
+      }                                                                      \
+    }                                                                        \
+    template <typename T>                                                    \
     bool operator()(const T* const x, T* residual) const {
 // clang-format off
 
@@ -223,7 +224,7 @@
   const T x1 = x[0];
   const T x2 = x[1];
   const T x3 = x[2];
-  const T theta = (0.5 / M_PI)  * atan(x2 / x1) + (x1 > 0.0 ? 0.0 : 0.5);
+  const T theta = (0.5 / constants::pi)  * atan(x2 / x1) + (x1 > 0.0 ? 0.0 : 0.5);
   residual[0] = 10.0 * (x3 - 10.0 * theta);
   residual[1] = 10.0 * (sqrt(x1 * x1 + x2 * x2) - 1.0);
   residual[2] = x3;
@@ -548,7 +549,7 @@
   }
 
   Problem problem;
-  problem.AddResidualBlock(TestProblem::Create(), NULL, x);
+  problem.AddResidualBlock(TestProblem::Create(), nullptr, x);
   double optimal_cost = TestProblem::unconstrained_optimal_cost;
 
   if (is_constrained) {
@@ -580,13 +581,11 @@
   return success;
 }
 
-}  // namespace examples
-}  // namespace ceres
+}  // namespace ceres::examples
 
 int main(int argc, char** argv) {
   GFLAGS_NAMESPACE::ParseCommandLineFlags(&argc, &argv, true);
   google::InitGoogleLogging(argv[0]);
-
   using ceres::examples::Solve;
 
   int unconstrained_problems = 0;
@@ -597,7 +596,8 @@
 
 #define UNCONSTRAINED_SOLVE(n)                              \
   ss << "Unconstrained Problem " << n << " : ";             \
-  if (FLAGS_problem == #n || FLAGS_problem == "all") {      \
+  if (CERES_GET_FLAG(FLAGS_problem) == #n ||                \
+      CERES_GET_FLAG(FLAGS_problem) == "all") {             \
     unconstrained_problems += 3;                            \
     if (Solve<ceres::examples::TestProblem##n>(false, 0)) { \
       unconstrained_successes += 1;                         \
@@ -645,7 +645,8 @@
 
 #define CONSTRAINED_SOLVE(n)                               \
   ss << "Constrained Problem " << n << " : ";              \
-  if (FLAGS_problem == #n || FLAGS_problem == "all") {     \
+  if (CERES_GET_FLAG(FLAGS_problem) == #n ||               \
+      CERES_GET_FLAG(FLAGS_problem) == "all") {            \
     constrained_problems += 1;                             \
     if (Solve<ceres::examples::TestProblem##n>(true, 0)) { \
       constrained_successes += 1;                          \
diff --git a/examples/nist.cc b/examples/nist.cc
index 977b69d..b92c918 100644
--- a/examples/nist.cc
+++ b/examples/nist.cc
@@ -1,5 +1,5 @@
 // Ceres Solver - A fast non-linear least squares minimizer
-// Copyright 2017 Google Inc. All rights reserved.
+// Copyright 2023 Google Inc. All rights reserved.
 // http://ceres-solver.org/
 //
 // Redistribution and use in source and binary forms, with or without
@@ -71,9 +71,12 @@
 // Average LRE     2.3      4.3       4.0  6.8      4.4    9.4
 //      Winner       0        0         5   11        2     41
 
+#include <cstdlib>
 #include <fstream>
 #include <iostream>
 #include <iterator>
+#include <string>
+#include <vector>
 
 #include "Eigen/Core"
 #include "ceres/ceres.h"
@@ -99,6 +102,9 @@
               "dense_qr",
               "Options are: sparse_cholesky, dense_qr, dense_normal_cholesky "
               "and cgnr");
+DEFINE_string(dense_linear_algebra_library,
+              "eigen",
+              "Options are: eigen, lapack, and cuda.");
 DEFINE_string(preconditioner, "jacobi", "Options are: identity, jacobi");
 DEFINE_string(line_search,
               "wolfe",
@@ -114,7 +120,7 @@
              "Maximum number of restarts of line search direction algorithm.");
 DEFINE_string(line_search_interpolation,
               "cubic",
-              "Degree of polynomial aproximation in line search, choices are: "
+              "Degree of polynomial approximation in line search, choices are: "
               "bisection, quadratic & cubic.");
 DEFINE_int32(lbfgs_rank,
              20,
@@ -148,26 +154,18 @@
              3,
              "Maximal number of Ridders extrapolations.");
 
-namespace ceres {
-namespace examples {
+namespace ceres::examples {
 namespace {
 
 using Eigen::Dynamic;
 using Eigen::RowMajor;
-typedef Eigen::Matrix<double, Dynamic, 1> Vector;
-typedef Eigen::Matrix<double, Dynamic, Dynamic, RowMajor> Matrix;
+using Vector = Eigen::Matrix<double, Dynamic, 1>;
+using Matrix = Eigen::Matrix<double, Dynamic, Dynamic, RowMajor>;
 
-using std::atof;
-using std::atoi;
-using std::cout;
-using std::ifstream;
-using std::string;
-using std::vector;
-
-void SplitStringUsingChar(const string& full,
+void SplitStringUsingChar(const std::string& full,
                           const char delim,
-                          vector<string>* result) {
-  std::back_insert_iterator<vector<string>> it(*result);
+                          std::vector<std::string>* result) {
+  std::back_insert_iterator<std::vector<std::string>> it(*result);
 
   const char* p = full.data();
   const char* end = p + full.size();
@@ -177,22 +175,22 @@
     } else {
       const char* start = p;
       while (++p != end && *p != delim) {
-        // Skip to the next occurence of the delimiter.
+        // Skip to the next occurrence of the delimiter.
       }
-      *it++ = string(start, p - start);
+      *it++ = std::string(start, p - start);
     }
   }
 }
 
-bool GetAndSplitLine(ifstream& ifs, vector<string>* pieces) {
+bool GetAndSplitLine(std::ifstream& ifs, std::vector<std::string>* pieces) {
   pieces->clear();
   char buf[256];
   ifs.getline(buf, 256);
-  SplitStringUsingChar(string(buf), ' ', pieces);
+  SplitStringUsingChar(std::string(buf), ' ', pieces);
   return true;
 }
 
-void SkipLines(ifstream& ifs, int num_lines) {
+void SkipLines(std::ifstream& ifs, int num_lines) {
   char buf[256];
   for (int i = 0; i < num_lines; ++i) {
     ifs.getline(buf, 256);
@@ -201,24 +199,24 @@
 
 class NISTProblem {
  public:
-  explicit NISTProblem(const string& filename) {
-    ifstream ifs(filename.c_str(), ifstream::in);
+  explicit NISTProblem(const std::string& filename) {
+    std::ifstream ifs(filename.c_str(), std::ifstream::in);
     CHECK(ifs) << "Unable to open : " << filename;
 
-    vector<string> pieces;
+    std::vector<std::string> pieces;
     SkipLines(ifs, 24);
     GetAndSplitLine(ifs, &pieces);
-    const int kNumResponses = atoi(pieces[1].c_str());
+    const int kNumResponses = std::atoi(pieces[1].c_str());
 
     GetAndSplitLine(ifs, &pieces);
-    const int kNumPredictors = atoi(pieces[0].c_str());
+    const int kNumPredictors = std::atoi(pieces[0].c_str());
 
     GetAndSplitLine(ifs, &pieces);
-    const int kNumObservations = atoi(pieces[0].c_str());
+    const int kNumObservations = std::atoi(pieces[0].c_str());
 
     SkipLines(ifs, 4);
     GetAndSplitLine(ifs, &pieces);
-    const int kNumParameters = atoi(pieces[0].c_str());
+    const int kNumParameters = std::atoi(pieces[0].c_str());
     SkipLines(ifs, 8);
 
     // Get the first line of initial and final parameter values to
@@ -234,24 +232,26 @@
     // Parse the line for parameter b1.
     int parameter_id = 0;
     for (int i = 0; i < kNumTries; ++i) {
-      initial_parameters_(i, parameter_id) = atof(pieces[i + 2].c_str());
+      initial_parameters_(i, parameter_id) = std::atof(pieces[i + 2].c_str());
     }
-    final_parameters_(0, parameter_id) = atof(pieces[2 + kNumTries].c_str());
+    final_parameters_(0, parameter_id) =
+        std::atof(pieces[2 + kNumTries].c_str());
 
     // Parse the remaining parameter lines.
     for (int parameter_id = 1; parameter_id < kNumParameters; ++parameter_id) {
       GetAndSplitLine(ifs, &pieces);
       // b2, b3, ....
       for (int i = 0; i < kNumTries; ++i) {
-        initial_parameters_(i, parameter_id) = atof(pieces[i + 2].c_str());
+        initial_parameters_(i, parameter_id) = std::atof(pieces[i + 2].c_str());
       }
-      final_parameters_(0, parameter_id) = atof(pieces[2 + kNumTries].c_str());
+      final_parameters_(0, parameter_id) =
+          std::atof(pieces[2 + kNumTries].c_str());
     }
 
-    // Certfied cost
+    // Certified cost
     SkipLines(ifs, 1);
     GetAndSplitLine(ifs, &pieces);
-    certified_cost_ = atof(pieces[4].c_str()) / 2.0;
+    certified_cost_ = std::atof(pieces[4].c_str()) / 2.0;
 
     // Read the observations.
     SkipLines(ifs, 18 - kNumParameters);
@@ -259,12 +259,12 @@
       GetAndSplitLine(ifs, &pieces);
       // Response.
       for (int j = 0; j < kNumResponses; ++j) {
-        response_(i, j) = atof(pieces[j].c_str());
+        response_(i, j) = std::atof(pieces[j].c_str());
       }
 
       // Predictor variables.
       for (int j = 0; j < kNumPredictors; ++j) {
-        predictor_(i, j) = atof(pieces[j + kNumResponses].c_str());
+        predictor_(i, j) = std::atof(pieces[j + kNumResponses].c_str());
       }
     }
   }
@@ -455,47 +455,57 @@
 // clang-format on
 
 static void SetNumericDiffOptions(ceres::NumericDiffOptions* options) {
-  options->max_num_ridders_extrapolations = FLAGS_ridders_extrapolations;
-  options->ridders_relative_initial_step_size = FLAGS_ridders_step_size;
+  options->max_num_ridders_extrapolations =
+      CERES_GET_FLAG(FLAGS_ridders_extrapolations);
+  options->ridders_relative_initial_step_size =
+      CERES_GET_FLAG(FLAGS_ridders_step_size);
 }
 
 void SetMinimizerOptions(ceres::Solver::Options* options) {
-  CHECK(
-      ceres::StringToMinimizerType(FLAGS_minimizer, &options->minimizer_type));
-  CHECK(ceres::StringToLinearSolverType(FLAGS_linear_solver,
+  CHECK(ceres::StringToMinimizerType(CERES_GET_FLAG(FLAGS_minimizer),
+                                     &options->minimizer_type));
+  CHECK(ceres::StringToLinearSolverType(CERES_GET_FLAG(FLAGS_linear_solver),
                                         &options->linear_solver_type));
-  CHECK(ceres::StringToPreconditionerType(FLAGS_preconditioner,
+  CHECK(StringToDenseLinearAlgebraLibraryType(
+      CERES_GET_FLAG(FLAGS_dense_linear_algebra_library),
+      &options->dense_linear_algebra_library_type));
+  CHECK(ceres::StringToPreconditionerType(CERES_GET_FLAG(FLAGS_preconditioner),
                                           &options->preconditioner_type));
   CHECK(ceres::StringToTrustRegionStrategyType(
-      FLAGS_trust_region_strategy, &options->trust_region_strategy_type));
-  CHECK(ceres::StringToDoglegType(FLAGS_dogleg, &options->dogleg_type));
+      CERES_GET_FLAG(FLAGS_trust_region_strategy),
+      &options->trust_region_strategy_type));
+  CHECK(ceres::StringToDoglegType(CERES_GET_FLAG(FLAGS_dogleg),
+                                  &options->dogleg_type));
   CHECK(ceres::StringToLineSearchDirectionType(
-      FLAGS_line_search_direction, &options->line_search_direction_type));
-  CHECK(ceres::StringToLineSearchType(FLAGS_line_search,
+      CERES_GET_FLAG(FLAGS_line_search_direction),
+      &options->line_search_direction_type));
+  CHECK(ceres::StringToLineSearchType(CERES_GET_FLAG(FLAGS_line_search),
                                       &options->line_search_type));
   CHECK(ceres::StringToLineSearchInterpolationType(
-      FLAGS_line_search_interpolation,
+      CERES_GET_FLAG(FLAGS_line_search_interpolation),
       &options->line_search_interpolation_type));
 
-  options->max_num_iterations = FLAGS_num_iterations;
-  options->use_nonmonotonic_steps = FLAGS_nonmonotonic_steps;
-  options->initial_trust_region_radius = FLAGS_initial_trust_region_radius;
-  options->max_lbfgs_rank = FLAGS_lbfgs_rank;
-  options->line_search_sufficient_function_decrease = FLAGS_sufficient_decrease;
+  options->max_num_iterations = CERES_GET_FLAG(FLAGS_num_iterations);
+  options->use_nonmonotonic_steps = CERES_GET_FLAG(FLAGS_nonmonotonic_steps);
+  options->initial_trust_region_radius =
+      CERES_GET_FLAG(FLAGS_initial_trust_region_radius);
+  options->max_lbfgs_rank = CERES_GET_FLAG(FLAGS_lbfgs_rank);
+  options->line_search_sufficient_function_decrease =
+      CERES_GET_FLAG(FLAGS_sufficient_decrease);
   options->line_search_sufficient_curvature_decrease =
-      FLAGS_sufficient_curvature_decrease;
+      CERES_GET_FLAG(FLAGS_sufficient_curvature_decrease);
   options->max_num_line_search_step_size_iterations =
-      FLAGS_max_line_search_iterations;
+      CERES_GET_FLAG(FLAGS_max_line_search_iterations);
   options->max_num_line_search_direction_restarts =
-      FLAGS_max_line_search_restarts;
+      CERES_GET_FLAG(FLAGS_max_line_search_restarts);
   options->use_approximate_eigenvalue_bfgs_scaling =
-      FLAGS_approximate_eigenvalue_bfgs_scaling;
+      CERES_GET_FLAG(FLAGS_approximate_eigenvalue_bfgs_scaling);
   options->function_tolerance = std::numeric_limits<double>::epsilon();
   options->gradient_tolerance = std::numeric_limits<double>::epsilon();
   options->parameter_tolerance = std::numeric_limits<double>::epsilon();
 }
 
-string JoinPath(const string& dirname, const string& basename) {
+std::string JoinPath(const std::string& dirname, const std::string& basename) {
 #ifdef _WIN32
   static const char separator = '\\';
 #else
@@ -507,7 +517,7 @@
   } else if (dirname[dirname.size() - 1] == separator) {
     return dirname + basename;
   } else {
-    return dirname + string(&separator, 1) + basename;
+    return dirname + std::string(&separator, 1) + basename;
   }
 }
 
@@ -515,24 +525,24 @@
 CostFunction* CreateCostFunction(const Matrix& predictor,
                                  const Matrix& response,
                                  const int num_observations) {
-  Model* model = new Model(predictor.data(), response.data(), num_observations);
-  ceres::CostFunction* cost_function = NULL;
-  if (FLAGS_use_numeric_diff) {
+  auto* model = new Model(predictor.data(), response.data(), num_observations);
+  ceres::CostFunction* cost_function = nullptr;
+  if (CERES_GET_FLAG(FLAGS_use_numeric_diff)) {
     ceres::NumericDiffOptions options;
     SetNumericDiffOptions(&options);
-    if (FLAGS_numeric_diff_method == "central") {
+    if (CERES_GET_FLAG(FLAGS_numeric_diff_method) == "central") {
       cost_function = new NumericDiffCostFunction<Model,
                                                   ceres::CENTRAL,
                                                   ceres::DYNAMIC,
                                                   num_parameters>(
           model, ceres::TAKE_OWNERSHIP, num_observations, options);
-    } else if (FLAGS_numeric_diff_method == "forward") {
+    } else if (CERES_GET_FLAG(FLAGS_numeric_diff_method) == "forward") {
       cost_function = new NumericDiffCostFunction<Model,
                                                   ceres::FORWARD,
                                                   ceres::DYNAMIC,
                                                   num_parameters>(
           model, ceres::TAKE_OWNERSHIP, num_observations, options);
-    } else if (FLAGS_numeric_diff_method == "ridders") {
+    } else if (CERES_GET_FLAG(FLAGS_numeric_diff_method) == "ridders") {
       cost_function = new NumericDiffCostFunction<Model,
                                                   ceres::RIDDERS,
                                                   ceres::DYNAMIC,
@@ -540,7 +550,7 @@
           model, ceres::TAKE_OWNERSHIP, num_observations, options);
     } else {
       LOG(ERROR) << "Invalid numeric diff method specified";
-      return 0;
+      return nullptr;
     }
   } else {
     cost_function =
@@ -571,8 +581,9 @@
 }
 
 template <typename Model, int num_parameters>
-int RegressionDriver(const string& filename) {
-  NISTProblem nist_problem(JoinPath(FLAGS_nist_data_dir, filename));
+int RegressionDriver(const std::string& filename) {
+  NISTProblem nist_problem(
+      JoinPath(CERES_GET_FLAG(FLAGS_nist_data_dir), filename));
   CHECK_EQ(num_parameters, nist_problem.num_parameters());
 
   Matrix predictor = nist_problem.predictor();
@@ -593,9 +604,10 @@
     double initial_cost;
     double final_cost;
 
-    if (!FLAGS_use_tiny_solver) {
+    if (!CERES_GET_FLAG(FLAGS_use_tiny_solver)) {
       ceres::Problem problem;
-      problem.AddResidualBlock(cost_function, NULL, initial_parameters.data());
+      problem.AddResidualBlock(
+          cost_function, nullptr, initial_parameters.data());
       ceres::Solver::Summary summary;
       ceres::Solver::Options options;
       SetMinimizerOptions(&options);
@@ -605,15 +617,15 @@
     } else {
       ceres::TinySolverCostFunctionAdapter<Eigen::Dynamic, num_parameters> cfa(
           *cost_function);
-      typedef ceres::TinySolver<
-          ceres::TinySolverCostFunctionAdapter<Eigen::Dynamic, num_parameters>>
-          Solver;
+      using Solver = ceres::TinySolver<
+          ceres::TinySolverCostFunctionAdapter<Eigen::Dynamic, num_parameters>>;
       Solver solver;
-      solver.options.max_num_iterations = FLAGS_num_iterations;
+      solver.options.max_num_iterations = CERES_GET_FLAG(FLAGS_num_iterations);
       solver.options.gradient_tolerance =
           std::numeric_limits<double>::epsilon();
       solver.options.parameter_tolerance =
           std::numeric_limits<double>::epsilon();
+      solver.options.function_tolerance = 0.0;
 
       Eigen::Matrix<double, num_parameters, 1> x;
       x = initial_parameters.transpose();
@@ -645,11 +657,11 @@
 }
 
 void SolveNISTProblems() {
-  if (FLAGS_nist_data_dir.empty()) {
+  if (CERES_GET_FLAG(FLAGS_nist_data_dir).empty()) {
     LOG(FATAL) << "Must specify the directory containing the NIST problems";
   }
 
-  cout << "Lower Difficulty\n";
+  std::cout << "Lower Difficulty\n";
   int easy_success = 0;
   easy_success += RegressionDriver<Misra1a, 2>("Misra1a.dat");
   easy_success += RegressionDriver<Chwirut, 3>("Chwirut1.dat");
@@ -660,7 +672,7 @@
   easy_success += RegressionDriver<DanWood, 2>("DanWood.dat");
   easy_success += RegressionDriver<Misra1b, 2>("Misra1b.dat");
 
-  cout << "\nMedium Difficulty\n";
+  std::cout << "\nMedium Difficulty\n";
   int medium_success = 0;
   medium_success += RegressionDriver<Kirby2, 5>("Kirby2.dat");
   medium_success += RegressionDriver<Hahn1, 7>("Hahn1.dat");
@@ -674,7 +686,7 @@
   medium_success += RegressionDriver<Roszman1, 4>("Roszman1.dat");
   medium_success += RegressionDriver<ENSO, 9>("ENSO.dat");
 
-  cout << "\nHigher Difficulty\n";
+  std::cout << "\nHigher Difficulty\n";
   int hard_success = 0;
   hard_success += RegressionDriver<MGH09, 4>("MGH09.dat");
   hard_success += RegressionDriver<Thurber, 7>("Thurber.dat");
@@ -685,17 +697,16 @@
   hard_success += RegressionDriver<Rat43, 4>("Rat43.dat");
   hard_success += RegressionDriver<Bennet5, 3>("Bennett5.dat");
 
-  cout << "\n";
-  cout << "Easy    : " << easy_success << "/16\n";
-  cout << "Medium  : " << medium_success << "/22\n";
-  cout << "Hard    : " << hard_success << "/16\n";
-  cout << "Total   : " << easy_success + medium_success + hard_success
-       << "/54\n";
+  std::cout << "\n";
+  std::cout << "Easy    : " << easy_success << "/16\n";
+  std::cout << "Medium  : " << medium_success << "/22\n";
+  std::cout << "Hard    : " << hard_success << "/16\n";
+  std::cout << "Total   : " << easy_success + medium_success + hard_success
+            << "/54\n";
 }
 
 }  // namespace
-}  // namespace examples
-}  // namespace ceres
+}  // namespace ceres::examples
 
 int main(int argc, char** argv) {
   GFLAGS_NAMESPACE::ParseCommandLineFlags(&argc, &argv, true);
diff --git a/examples/pgm_image.h b/examples/pgm_image.h
index 3d2df63..033ab4d 100644
--- a/examples/pgm_image.h
+++ b/examples/pgm_image.h
@@ -1,5 +1,5 @@
 // Ceres Solver - A fast non-linear least squares minimizer
-// Copyright 2015 Google Inc. All rights reserved.
+// Copyright 2023 Google Inc. All rights reserved.
 // http://ceres-solver.org/
 //
 // Redistribution and use in source and binary forms, with or without
@@ -43,8 +43,7 @@
 
 #include "glog/logging.h"
 
-namespace ceres {
-namespace examples {
+namespace ceres::examples {
 
 template <typename Real>
 class PGMImage {
@@ -311,7 +310,6 @@
   return data_;
 }
 
-}  // namespace examples
-}  // namespace ceres
+}  // namespace ceres::examples
 
 #endif  // CERES_EXAMPLES_PGM_IMAGE_H_
diff --git a/examples/powell.cc b/examples/powell.cc
index c75ad24..a4ca1b7 100644
--- a/examples/powell.cc
+++ b/examples/powell.cc
@@ -1,5 +1,5 @@
 // Ceres Solver - A fast non-linear least squares minimizer
-// Copyright 2015 Google Inc. All rights reserved.
+// Copyright 2023 Google Inc. All rights reserved.
 // http://ceres-solver.org/
 //
 // Redistribution and use in source and binary forms, with or without
@@ -50,12 +50,6 @@
 #include "gflags/gflags.h"
 #include "glog/logging.h"
 
-using ceres::AutoDiffCostFunction;
-using ceres::CostFunction;
-using ceres::Problem;
-using ceres::Solve;
-using ceres::Solver;
-
 struct F1 {
   template <typename T>
   bool operator()(const T* const x1, const T* const x2, T* residual) const {
@@ -105,24 +99,24 @@
   double x3 = 0.0;
   double x4 = 1.0;
 
-  Problem problem;
-  // Add residual terms to the problem using the using the autodiff
+  ceres::Problem problem;
+  // Add residual terms to the problem using the autodiff
   // wrapper to get the derivatives automatically. The parameters, x1 through
   // x4, are modified in place.
   problem.AddResidualBlock(
-      new AutoDiffCostFunction<F1, 1, 1, 1>(new F1), NULL, &x1, &x2);
+      new ceres::AutoDiffCostFunction<F1, 1, 1, 1>(), nullptr, &x1, &x2);
   problem.AddResidualBlock(
-      new AutoDiffCostFunction<F2, 1, 1, 1>(new F2), NULL, &x3, &x4);
+      new ceres::AutoDiffCostFunction<F2, 1, 1, 1>(), nullptr, &x3, &x4);
   problem.AddResidualBlock(
-      new AutoDiffCostFunction<F3, 1, 1, 1>(new F3), NULL, &x2, &x3);
+      new ceres::AutoDiffCostFunction<F3, 1, 1, 1>(), nullptr, &x2, &x3);
   problem.AddResidualBlock(
-      new AutoDiffCostFunction<F4, 1, 1, 1>(new F4), NULL, &x1, &x4);
+      new ceres::AutoDiffCostFunction<F4, 1, 1, 1>(), nullptr, &x1, &x4);
 
-  Solver::Options options;
-  LOG_IF(
-      FATAL,
-      !ceres::StringToMinimizerType(FLAGS_minimizer, &options.minimizer_type))
-      << "Invalid minimizer: " << FLAGS_minimizer
+  ceres::Solver::Options options;
+  LOG_IF(FATAL,
+         !ceres::StringToMinimizerType(CERES_GET_FLAG(FLAGS_minimizer),
+                                       &options.minimizer_type))
+      << "Invalid minimizer: " << CERES_GET_FLAG(FLAGS_minimizer)
       << ", valid options are: trust_region and line_search.";
 
   options.max_num_iterations = 100;
@@ -138,8 +132,8 @@
   // clang-format on
 
   // Run the solver!
-  Solver::Summary summary;
-  Solve(options, &problem, &summary);
+  ceres::Solver::Summary summary;
+  ceres::Solve(options, &problem, &summary);
 
   std::cout << summary.FullReport() << "\n";
   // clang-format off
diff --git a/examples/random.h b/examples/random.h
deleted file mode 100644
index ace0711..0000000
--- a/examples/random.h
+++ /dev/null
@@ -1,64 +0,0 @@
-// Ceres Solver - A fast non-linear least squares minimizer
-// Copyright 2015 Google Inc. All rights reserved.
-// http://ceres-solver.org/
-//
-// Redistribution and use in source and binary forms, with or without
-// modification, are permitted provided that the following conditions are met:
-//
-// * Redistributions of source code must retain the above copyright notice,
-//   this list of conditions and the following disclaimer.
-// * Redistributions in binary form must reproduce the above copyright notice,
-//   this list of conditions and the following disclaimer in the documentation
-//   and/or other materials provided with the distribution.
-// * Neither the name of Google Inc. nor the names of its contributors may be
-//   used to endorse or promote products derived from this software without
-//   specific prior written permission.
-//
-// THIS SOFTWARE IS PROVIDED BY THE COPYRIGHT HOLDERS AND CONTRIBUTORS "AS IS"
-// AND ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED TO, THE
-// IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR PURPOSE
-// ARE DISCLAIMED. IN NO EVENT SHALL THE COPYRIGHT OWNER OR CONTRIBUTORS BE
-// LIABLE FOR ANY DIRECT, INDIRECT, INCIDENTAL, SPECIAL, EXEMPLARY, OR
-// CONSEQUENTIAL DAMAGES (INCLUDING, BUT NOT LIMITED TO, PROCUREMENT OF
-// SUBSTITUTE GOODS OR SERVICES; LOSS OF USE, DATA, OR PROFITS; OR BUSINESS
-// INTERRUPTION) HOWEVER CAUSED AND ON ANY THEORY OF LIABILITY, WHETHER IN
-// CONTRACT, STRICT LIABILITY, OR TORT (INCLUDING NEGLIGENCE OR OTHERWISE)
-// ARISING IN ANY WAY OUT OF THE USE OF THIS SOFTWARE, EVEN IF ADVISED OF THE
-// POSSIBILITY OF SUCH DAMAGE.
-//
-// Author: sameeragarwal@google.com (Sameer Agarwal)
-
-#ifndef CERES_EXAMPLES_RANDOM_H_
-#define CERES_EXAMPLES_RANDOM_H_
-
-#include <math.h>
-#include <stdlib.h>
-
-namespace ceres {
-namespace examples {
-
-// Return a random number sampled from a uniform distribution in the range
-// [0,1].
-inline double RandDouble() {
-  double r = static_cast<double>(rand());
-  return r / RAND_MAX;
-}
-
-// Marsaglia Polar method for generation standard normal (pseudo)
-// random numbers http://en.wikipedia.org/wiki/Marsaglia_polar_method
-inline double RandNormal() {
-  double x1, x2, w;
-  do {
-    x1 = 2.0 * RandDouble() - 1.0;
-    x2 = 2.0 * RandDouble() - 1.0;
-    w = x1 * x1 + x2 * x2;
-  } while (w >= 1.0 || w == 0.0);
-
-  w = sqrt((-2.0 * log(w)) / w);
-  return x1 * w;
-}
-
-}  // namespace examples
-}  // namespace ceres
-
-#endif  // CERES_EXAMPLES_RANDOM_H_
diff --git a/examples/robot_pose_mle.cc b/examples/robot_pose_mle.cc
index ab9a098..cc60e14 100644
--- a/examples/robot_pose_mle.cc
+++ b/examples/robot_pose_mle.cc
@@ -1,5 +1,5 @@
 // Ceres Solver - A fast non-linear least squares minimizer
-// Copyright 2015 Google Inc. All rights reserved.
+// Copyright 2023 Google Inc. All rights reserved.
 // http://ceres-solver.org/
 //
 // Redistribution and use in source and binary forms, with or without
@@ -50,7 +50,7 @@
 //
 // There are two types of residuals in this problem:
 // 1) The OdometryConstraint residual, that accounts for the odometry readings
-//    between successive pose estimatess of the robot.
+//    between successive pose estimates of the robot.
 // 2) The RangeConstraint residual, that accounts for the errors in the observed
 //    range readings from each pose.
 //
@@ -97,14 +97,14 @@
 // timesteps 0 to i for that variable, both inclusive.
 //
 // Bayes' rule is used to derive eq. 3 from 2, and the independence of
-// odometry observations and range readings is expolited to derive 4 from 3.
+// odometry observations and range readings is exploited to derive 4 from 3.
 //
 // Thus, the Belief, up to scale, is factored as a product of a number of
 // terms, two for each pose, where for each pose term there is one term for the
 // range reading, P(y_i | u*_(0:i) and one term for the odometry reading,
 // P(u*_i | u_i) . Note that the term for the range reading is dependent on all
 // odometry values u*_(0:i), while the odometry term, P(u*_i | u_i) depends only
-// on a single value, u_i. Both the range reading as well as odoemtry
+// on a single value, u_i. Both the range reading as well as odometry
 // probability terms are modeled as the Normal distribution, and have the form:
 //
 // p(x) \propto \exp{-((x - x_mean) / x_stddev)^2}
@@ -123,30 +123,18 @@
 // variable, and will be computed by an AutoDiffCostFunction, while the term
 // for the range reading will depend on all previous odometry observations, and
 // will be computed by a DynamicAutoDiffCostFunction since the number of
-// odoemtry observations will only be known at run time.
+// odometry observations will only be known at run time.
 
-#include <math.h>
-
+#include <algorithm>
+#include <cmath>
 #include <cstdio>
+#include <random>
 #include <vector>
 
 #include "ceres/ceres.h"
 #include "ceres/dynamic_autodiff_cost_function.h"
 #include "gflags/gflags.h"
 #include "glog/logging.h"
-#include "random.h"
-
-using ceres::AutoDiffCostFunction;
-using ceres::CauchyLoss;
-using ceres::CostFunction;
-using ceres::DynamicAutoDiffCostFunction;
-using ceres::LossFunction;
-using ceres::Problem;
-using ceres::Solve;
-using ceres::Solver;
-using ceres::examples::RandNormal;
-using std::min;
-using std::vector;
 
 DEFINE_double(corridor_length,
               30.0,
@@ -169,7 +157,8 @@
 static constexpr int kStride = 10;
 
 struct OdometryConstraint {
-  typedef AutoDiffCostFunction<OdometryConstraint, 1, 1> OdometryCostFunction;
+  using OdometryCostFunction =
+      ceres::AutoDiffCostFunction<OdometryConstraint, 1, 1>;
 
   OdometryConstraint(double odometry_mean, double odometry_stddev)
       : odometry_mean(odometry_mean), odometry_stddev(odometry_stddev) {}
@@ -181,8 +170,8 @@
   }
 
   static OdometryCostFunction* Create(const double odometry_value) {
-    return new OdometryCostFunction(
-        new OdometryConstraint(odometry_value, FLAGS_odometry_stddev));
+    return new OdometryCostFunction(new OdometryConstraint(
+        odometry_value, CERES_GET_FLAG(FLAGS_odometry_stddev)));
   }
 
   const double odometry_mean;
@@ -190,8 +179,8 @@
 };
 
 struct RangeConstraint {
-  typedef DynamicAutoDiffCostFunction<RangeConstraint, kStride>
-      RangeCostFunction;
+  using RangeCostFunction =
+      ceres::DynamicAutoDiffCostFunction<RangeConstraint, kStride>;
 
   RangeConstraint(int pose_index,
                   double range_reading,
@@ -217,11 +206,14 @@
   // conveniently add to a ceres problem.
   static RangeCostFunction* Create(const int pose_index,
                                    const double range_reading,
-                                   vector<double>* odometry_values,
-                                   vector<double*>* parameter_blocks) {
-    RangeConstraint* constraint = new RangeConstraint(
-        pose_index, range_reading, FLAGS_range_stddev, FLAGS_corridor_length);
-    RangeCostFunction* cost_function = new RangeCostFunction(constraint);
+                                   std::vector<double>* odometry_values,
+                                   std::vector<double*>* parameter_blocks) {
+    auto* constraint =
+        new RangeConstraint(pose_index,
+                            range_reading,
+                            CERES_GET_FLAG(FLAGS_range_stddev),
+                            CERES_GET_FLAG(FLAGS_corridor_length));
+    auto* cost_function = new RangeCostFunction(constraint);
     // Add all the parameter blocks that affect this constraint.
     parameter_blocks->clear();
     for (int i = 0; i <= pose_index; ++i) {
@@ -240,37 +232,45 @@
 
 namespace {
 
-void SimulateRobot(vector<double>* odometry_values,
-                   vector<double>* range_readings) {
+void SimulateRobot(std::vector<double>* odometry_values,
+                   std::vector<double>* range_readings) {
   const int num_steps =
-      static_cast<int>(ceil(FLAGS_corridor_length / FLAGS_pose_separation));
+      static_cast<int>(ceil(CERES_GET_FLAG(FLAGS_corridor_length) /
+                            CERES_GET_FLAG(FLAGS_pose_separation)));
+  std::mt19937 prng;
+  std::normal_distribution<double> odometry_noise(
+      0.0, CERES_GET_FLAG(FLAGS_odometry_stddev));
+  std::normal_distribution<double> range_noise(
+      0.0, CERES_GET_FLAG(FLAGS_range_stddev));
 
   // The robot starts out at the origin.
   double robot_location = 0.0;
   for (int i = 0; i < num_steps; ++i) {
     const double actual_odometry_value =
-        min(FLAGS_pose_separation, FLAGS_corridor_length - robot_location);
+        std::min(CERES_GET_FLAG(FLAGS_pose_separation),
+                 CERES_GET_FLAG(FLAGS_corridor_length) - robot_location);
     robot_location += actual_odometry_value;
-    const double actual_range = FLAGS_corridor_length - robot_location;
+    const double actual_range =
+        CERES_GET_FLAG(FLAGS_corridor_length) - robot_location;
     const double observed_odometry =
-        RandNormal() * FLAGS_odometry_stddev + actual_odometry_value;
-    const double observed_range =
-        RandNormal() * FLAGS_range_stddev + actual_range;
+        actual_odometry_value + odometry_noise(prng);
+    const double observed_range = actual_range + range_noise(prng);
     odometry_values->push_back(observed_odometry);
     range_readings->push_back(observed_range);
   }
 }
 
-void PrintState(const vector<double>& odometry_readings,
-                const vector<double>& range_readings) {
+void PrintState(const std::vector<double>& odometry_readings,
+                const std::vector<double>& range_readings) {
   CHECK_EQ(odometry_readings.size(), range_readings.size());
   double robot_location = 0.0;
   printf("pose: location     odom    range  r.error  o.error\n");
   for (int i = 0; i < odometry_readings.size(); ++i) {
     robot_location += odometry_readings[i];
-    const double range_error =
-        robot_location + range_readings[i] - FLAGS_corridor_length;
-    const double odometry_error = FLAGS_pose_separation - odometry_readings[i];
+    const double range_error = robot_location + range_readings[i] -
+                               CERES_GET_FLAG(FLAGS_corridor_length);
+    const double odometry_error =
+        CERES_GET_FLAG(FLAGS_pose_separation) - odometry_readings[i];
     printf("%4d: %8.3f %8.3f %8.3f %8.3f %8.3f\n",
            static_cast<int>(i),
            robot_location,
@@ -287,13 +287,13 @@
   google::InitGoogleLogging(argv[0]);
   GFLAGS_NAMESPACE::ParseCommandLineFlags(&argc, &argv, true);
   // Make sure that the arguments parsed are all positive.
-  CHECK_GT(FLAGS_corridor_length, 0.0);
-  CHECK_GT(FLAGS_pose_separation, 0.0);
-  CHECK_GT(FLAGS_odometry_stddev, 0.0);
-  CHECK_GT(FLAGS_range_stddev, 0.0);
+  CHECK_GT(CERES_GET_FLAG(FLAGS_corridor_length), 0.0);
+  CHECK_GT(CERES_GET_FLAG(FLAGS_pose_separation), 0.0);
+  CHECK_GT(CERES_GET_FLAG(FLAGS_odometry_stddev), 0.0);
+  CHECK_GT(CERES_GET_FLAG(FLAGS_range_stddev), 0.0);
 
-  vector<double> odometry_values;
-  vector<double> range_readings;
+  std::vector<double> odometry_values;
+  std::vector<double> range_readings;
   SimulateRobot(&odometry_values, &range_readings);
 
   printf("Initial values:\n");
@@ -303,25 +303,25 @@
   for (int i = 0; i < odometry_values.size(); ++i) {
     // Create and add a DynamicAutoDiffCostFunction for the RangeConstraint from
     // pose i.
-    vector<double*> parameter_blocks;
+    std::vector<double*> parameter_blocks;
     RangeConstraint::RangeCostFunction* range_cost_function =
         RangeConstraint::Create(
             i, range_readings[i], &odometry_values, &parameter_blocks);
-    problem.AddResidualBlock(range_cost_function, NULL, parameter_blocks);
+    problem.AddResidualBlock(range_cost_function, nullptr, parameter_blocks);
 
     // Create and add an AutoDiffCostFunction for the OdometryConstraint for
     // pose i.
     problem.AddResidualBlock(OdometryConstraint::Create(odometry_values[i]),
-                             NULL,
+                             nullptr,
                              &(odometry_values[i]));
   }
 
   ceres::Solver::Options solver_options;
   solver_options.minimizer_progress_to_stdout = true;
 
-  Solver::Summary summary;
+  ceres::Solver::Summary summary;
   printf("Solving...\n");
-  Solve(solver_options, &problem, &summary);
+  ceres::Solve(solver_options, &problem, &summary);
   printf("Done.\n");
   std::cout << summary.FullReport() << "\n";
   printf("Final values:\n");
diff --git a/examples/robust_curve_fitting.cc b/examples/robust_curve_fitting.cc
index 9b526c5..e08b0df 100644
--- a/examples/robust_curve_fitting.cc
+++ b/examples/robust_curve_fitting.cc
@@ -1,5 +1,5 @@
 // Ceres Solver - A fast non-linear least squares minimizer
-// Copyright 2015 Google Inc. All rights reserved.
+// Copyright 2023 Google Inc. All rights reserved.
 // http://ceres-solver.org/
 //
 // Redistribution and use in source and binary forms, with or without
@@ -27,6 +27,12 @@
 // POSSIBILITY OF SUCH DAMAGE.
 //
 // Author: sameeragarwal@google.com (Sameer Agarwal)
+//
+// This example fits the curve f(x;m,c) = e^(m * x + c) to data. However unlike
+// the data in curve_fitting.cc, the data here has outliers in it, so minimizing
+// the sum squared loss will result in a bad fit. So this example illustrates
+// the use of a robust loss function (CauchyLoss) to reduce the influence of the
+// outliers on the fit.
 
 #include "ceres/ceres.h"
 #include "glog/logging.h"
@@ -115,13 +121,6 @@
 };
 // clang-format on
 
-using ceres::AutoDiffCostFunction;
-using ceres::CauchyLoss;
-using ceres::CostFunction;
-using ceres::Problem;
-using ceres::Solve;
-using ceres::Solver;
-
 struct ExponentialResidual {
   ExponentialResidual(double x, double y) : x_(x), y_(y) {}
 
@@ -139,25 +138,28 @@
 int main(int argc, char** argv) {
   google::InitGoogleLogging(argv[0]);
 
-  double m = 0.0;
-  double c = 0.0;
+  const double initial_m = 0.0;
+  const double initial_c = 0.0;
+  double m = initial_m;
+  double c = initial_c;
 
-  Problem problem;
+  ceres::Problem problem;
   for (int i = 0; i < kNumObservations; ++i) {
-    CostFunction* cost_function =
-        new AutoDiffCostFunction<ExponentialResidual, 1, 1, 1>(
-            new ExponentialResidual(data[2 * i], data[2 * i + 1]));
-    problem.AddResidualBlock(cost_function, new CauchyLoss(0.5), &m, &c);
+    ceres::CostFunction* cost_function =
+        new ceres::AutoDiffCostFunction<ExponentialResidual, 1, 1, 1>(
+            data[2 * i], data[2 * i + 1]);
+    problem.AddResidualBlock(cost_function, new ceres::CauchyLoss(0.5), &m, &c);
   }
 
-  Solver::Options options;
+  ceres::Solver::Options options;
+  options.max_num_iterations = 25;
   options.linear_solver_type = ceres::DENSE_QR;
   options.minimizer_progress_to_stdout = true;
 
-  Solver::Summary summary;
-  Solve(options, &problem, &summary);
+  ceres::Solver::Summary summary;
+  ceres::Solve(options, &problem, &summary);
   std::cout << summary.BriefReport() << "\n";
-  std::cout << "Initial m: " << 0.0 << " c: " << 0.0 << "\n";
+  std::cout << "Initial m: " << initial_m << " c: " << initial_c << "\n";
   std::cout << "Final   m: " << m << " c: " << c << "\n";
   return 0;
 }
diff --git a/examples/rosenbrock.cc b/examples/rosenbrock.cc
index 1b9aef6..a382ccd 100644
--- a/examples/rosenbrock.cc
+++ b/examples/rosenbrock.cc
@@ -1,5 +1,5 @@
 // Ceres Solver - A fast non-linear least squares minimizer
-// Copyright 2015 Google Inc. All rights reserved.
+// Copyright 2023 Google Inc. All rights reserved.
 // http://ceres-solver.org/
 //
 // Redistribution and use in source and binary forms, with or without
@@ -27,30 +27,29 @@
 // POSSIBILITY OF SUCH DAMAGE.
 //
 // Author: sameeragarwal@google.com (Sameer Agarwal)
+//
+// Example of minimizing the Rosenbrock function
+// (https://en.wikipedia.org/wiki/Rosenbrock_function) using
+// GradientProblemSolver using automatically computed derivatives.
 
 #include "ceres/ceres.h"
 #include "glog/logging.h"
 
 // f(x,y) = (1-x)^2 + 100(y - x^2)^2;
-class Rosenbrock : public ceres::FirstOrderFunction {
- public:
-  virtual ~Rosenbrock() {}
-
-  virtual bool Evaluate(const double* parameters,
-                        double* cost,
-                        double* gradient) const {
-    const double x = parameters[0];
-    const double y = parameters[1];
-
+struct Rosenbrock {
+  template <typename T>
+  bool operator()(const T* parameters, T* cost) const {
+    const T x = parameters[0];
+    const T y = parameters[1];
     cost[0] = (1.0 - x) * (1.0 - x) + 100.0 * (y - x * x) * (y - x * x);
-    if (gradient != NULL) {
-      gradient[0] = -2.0 * (1.0 - x) - 200.0 * (y - x * x) * 2.0 * x;
-      gradient[1] = 200.0 * (y - x * x);
-    }
     return true;
   }
 
-  virtual int NumParameters() const { return 2; }
+  static ceres::FirstOrderFunction* Create() {
+    constexpr int kNumParameters = 2;
+    return new ceres::AutoDiffFirstOrderFunction<Rosenbrock, kNumParameters>(
+        new Rosenbrock);
+  }
 };
 
 int main(int argc, char** argv) {
@@ -62,7 +61,7 @@
   options.minimizer_progress_to_stdout = true;
 
   ceres::GradientProblemSolver::Summary summary;
-  ceres::GradientProblem problem(new Rosenbrock());
+  ceres::GradientProblem problem(Rosenbrock::Create());
   ceres::Solve(options, problem, parameters, &summary);
 
   std::cout << summary.FullReport() << "\n";
diff --git a/examples/rosenbrock_analytic_diff.cc b/examples/rosenbrock_analytic_diff.cc
new file mode 100644
index 0000000..65e49eb
--- /dev/null
+++ b/examples/rosenbrock_analytic_diff.cc
@@ -0,0 +1,77 @@
+// Ceres Solver - A fast non-linear least squares minimizer
+// Copyright 2023 Google Inc. All rights reserved.
+// http://ceres-solver.org/
+//
+// Redistribution and use in source and binary forms, with or without
+// modification, are permitted provided that the following conditions are met:
+//
+// * Redistributions of source code must retain the above copyright notice,
+//   this list of conditions and the following disclaimer.
+// * Redistributions in binary form must reproduce the above copyright notice,
+//   this list of conditions and the following disclaimer in the documentation
+//   and/or other materials provided with the distribution.
+// * Neither the name of Google Inc. nor the names of its contributors may be
+//   used to endorse or promote products derived from this software without
+//   specific prior written permission.
+//
+// THIS SOFTWARE IS PROVIDED BY THE COPYRIGHT HOLDERS AND CONTRIBUTORS "AS IS"
+// AND ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED TO, THE
+// IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR PURPOSE
+// ARE DISCLAIMED. IN NO EVENT SHALL THE COPYRIGHT OWNER OR CONTRIBUTORS BE
+// LIABLE FOR ANY DIRECT, INDIRECT, INCIDENTAL, SPECIAL, EXEMPLARY, OR
+// CONSEQUENTIAL DAMAGES (INCLUDING, BUT NOT LIMITED TO, PROCUREMENT OF
+// SUBSTITUTE GOODS OR SERVICES; LOSS OF USE, DATA, OR PROFITS; OR BUSINESS
+// INTERRUPTION) HOWEVER CAUSED AND ON ANY THEORY OF LIABILITY, WHETHER IN
+// CONTRACT, STRICT LIABILITY, OR TORT (INCLUDING NEGLIGENCE OR OTHERWISE)
+// ARISING IN ANY WAY OUT OF THE USE OF THIS SOFTWARE, EVEN IF ADVISED OF THE
+// POSSIBILITY OF SUCH DAMAGE.
+//
+// Author: sameeragarwal@google.com (Sameer Agarwal)
+//
+// Example of minimizing the Rosenbrock function
+// (https://en.wikipedia.org/wiki/Rosenbrock_function) using
+// GradientProblemSolver using analytic derivatives.
+
+#include "ceres/ceres.h"
+#include "glog/logging.h"
+
+// f(x,y) = (1-x)^2 + 100(y - x^2)^2;
+class Rosenbrock final : public ceres::FirstOrderFunction {
+ public:
+  bool Evaluate(const double* parameters,
+                double* cost,
+                double* gradient) const override {
+    const double x = parameters[0];
+    const double y = parameters[1];
+
+    cost[0] = (1.0 - x) * (1.0 - x) + 100.0 * (y - x * x) * (y - x * x);
+
+    if (gradient) {
+      gradient[0] = -2.0 * (1.0 - x) - 200.0 * (y - x * x) * 2.0 * x;
+      gradient[1] = 200.0 * (y - x * x);
+    }
+
+    return true;
+  }
+
+  int NumParameters() const override { return 2; }
+};
+
+int main(int argc, char** argv) {
+  google::InitGoogleLogging(argv[0]);
+
+  double parameters[2] = {-1.2, 1.0};
+
+  ceres::GradientProblemSolver::Options options;
+  options.minimizer_progress_to_stdout = true;
+
+  ceres::GradientProblemSolver::Summary summary;
+  ceres::GradientProblem problem(new Rosenbrock());
+  ceres::Solve(options, problem, parameters, &summary);
+
+  std::cout << summary.FullReport() << "\n";
+  std::cout << "Initial x: " << -1.2 << " y: " << 1.0 << "\n";
+  std::cout << "Final   x: " << parameters[0] << " y: " << parameters[1]
+            << "\n";
+  return 0;
+}
diff --git a/examples/rosenbrock_numeric_diff.cc b/examples/rosenbrock_numeric_diff.cc
new file mode 100644
index 0000000..a711b2f
--- /dev/null
+++ b/examples/rosenbrock_numeric_diff.cc
@@ -0,0 +1,74 @@
+// Ceres Solver - A fast non-linear least squares minimizer
+// Copyright 2023 Google Inc. All rights reserved.
+// http://ceres-solver.org/
+//
+// Redistribution and use in source and binary forms, with or without
+// modification, are permitted provided that the following conditions are met:
+//
+// * Redistributions of source code must retain the above copyright notice,
+//   this list of conditions and the following disclaimer.
+// * Redistributions in binary form must reproduce the above copyright notice,
+//   this list of conditions and the following disclaimer in the documentation
+//   and/or other materials provided with the distribution.
+// * Neither the name of Google Inc. nor the names of its contributors may be
+//   used to endorse or promote products derived from this software without
+//   specific prior written permission.
+//
+// THIS SOFTWARE IS PROVIDED BY THE COPYRIGHT HOLDERS AND CONTRIBUTORS "AS IS"
+// AND ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED TO, THE
+// IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR PURPOSE
+// ARE DISCLAIMED. IN NO EVENT SHALL THE COPYRIGHT OWNER OR CONTRIBUTORS BE
+// LIABLE FOR ANY DIRECT, INDIRECT, INCIDENTAL, SPECIAL, EXEMPLARY, OR
+// CONSEQUENTIAL DAMAGES (INCLUDING, BUT NOT LIMITED TO, PROCUREMENT OF
+// SUBSTITUTE GOODS OR SERVICES; LOSS OF USE, DATA, OR PROFITS; OR BUSINESS
+// INTERRUPTION) HOWEVER CAUSED AND ON ANY THEORY OF LIABILITY, WHETHER IN
+// CONTRACT, STRICT LIABILITY, OR TORT (INCLUDING NEGLIGENCE OR OTHERWISE)
+// ARISING IN ANY WAY OUT OF THE USE OF THIS SOFTWARE, EVEN IF ADVISED OF THE
+// POSSIBILITY OF SUCH DAMAGE.
+//
+// Author: sameeragarwal@google.com (Sameer Agarwal)
+//
+// Example of minimizing the Rosenbrock function
+// (https://en.wikipedia.org/wiki/Rosenbrock_function) using
+// GradientProblemSolver using derivatives computed using numeric
+// differentiation.
+
+#include "ceres/ceres.h"
+#include "glog/logging.h"
+
+// f(x,y) = (1-x)^2 + 100(y - x^2)^2;
+struct Rosenbrock {
+  bool operator()(const double* parameters, double* cost) const {
+    const double x = parameters[0];
+    const double y = parameters[1];
+    cost[0] = (1.0 - x) * (1.0 - x) + 100.0 * (y - x * x) * (y - x * x);
+    return true;
+  }
+
+  static ceres::FirstOrderFunction* Create() {
+    constexpr int kNumParameters = 2;
+    return new ceres::NumericDiffFirstOrderFunction<Rosenbrock,
+                                                    ceres::CENTRAL,
+                                                    kNumParameters>(
+        new Rosenbrock);
+  }
+};
+
+int main(int argc, char** argv) {
+  google::InitGoogleLogging(argv[0]);
+
+  double parameters[2] = {-1.2, 1.0};
+
+  ceres::GradientProblemSolver::Options options;
+  options.minimizer_progress_to_stdout = true;
+
+  ceres::GradientProblemSolver::Summary summary;
+  ceres::GradientProblem problem(Rosenbrock::Create());
+  ceres::Solve(options, problem, parameters, &summary);
+
+  std::cout << summary.FullReport() << "\n";
+  std::cout << "Initial x: " << -1.2 << " y: " << 1.0 << "\n";
+  std::cout << "Final   x: " << parameters[0] << " y: " << parameters[1]
+            << "\n";
+  return 0;
+}
diff --git a/examples/sampled_function/CMakeLists.txt b/examples/sampled_function/CMakeLists.txt
index 8a17cad..1013753 100644
--- a/examples/sampled_function/CMakeLists.txt
+++ b/examples/sampled_function/CMakeLists.txt
@@ -1,5 +1,5 @@
 # Ceres Solver - A fast non-linear least squares minimizer
-# Copyright 2015 Google Inc. All rights reserved.
+# Copyright 2023 Google Inc. All rights reserved.
 # http://ceres-solver.org/
 #
 # Redistribution and use in source and binary forms, with or without
@@ -29,4 +29,4 @@
 # Author: vitus@google.com (Michael Vitus)
 
 add_executable(sampled_function sampled_function.cc)
-target_link_libraries(sampled_function Ceres::ceres)
+target_link_libraries(sampled_function PRIVATE Ceres::ceres)
diff --git a/examples/sampled_function/README.md b/examples/sampled_function/README.md
index ef1af43..5fde415 100644
--- a/examples/sampled_function/README.md
+++ b/examples/sampled_function/README.md
@@ -32,7 +32,7 @@
 
 ```c++
 bool Evaluate(double const* const* parameters, double* residuals, double** jacobians) const {
-  if (jacobians == NULL || jacobians[0] == NULL)
+  if (jacobians == nullptr || jacobians[0] == nullptr)
     interpolator_.Evaluate(parameters[0][0], residuals);
   else
     interpolator_.Evaluate(parameters[0][0], residuals, jacobians[0]);
diff --git a/examples/sampled_function/sampled_function.cc b/examples/sampled_function/sampled_function.cc
index e96018d..40e9c1f 100644
--- a/examples/sampled_function/sampled_function.cc
+++ b/examples/sampled_function/sampled_function.cc
@@ -1,5 +1,5 @@
 // Ceres Solver - A fast non-linear least squares minimizer
-// Copyright 2015 Google Inc. All rights reserved.
+// Copyright 2023 Google Inc. All rights reserved.
 // http://ceres-solver.org/
 //
 // Redistribution and use in source and binary forms, with or without
@@ -35,35 +35,27 @@
 #include "ceres/cubic_interpolation.h"
 #include "glog/logging.h"
 
-using ceres::AutoDiffCostFunction;
-using ceres::CostFunction;
-using ceres::CubicInterpolator;
-using ceres::Grid1D;
-using ceres::Problem;
-using ceres::Solve;
-using ceres::Solver;
+using Interpolator = ceres::CubicInterpolator<ceres::Grid1D<double>>;
 
 // A simple cost functor that interfaces an interpolated table of
 // values with automatic differentiation.
 struct InterpolatedCostFunctor {
-  explicit InterpolatedCostFunctor(
-      const CubicInterpolator<Grid1D<double>>& interpolator)
-      : interpolator_(interpolator) {}
+  explicit InterpolatedCostFunctor(const Interpolator& interpolator)
+      : interpolator(interpolator) {}
 
   template <typename T>
   bool operator()(const T* x, T* residuals) const {
-    interpolator_.Evaluate(*x, residuals);
+    interpolator.Evaluate(*x, residuals);
     return true;
   }
 
-  static CostFunction* Create(
-      const CubicInterpolator<Grid1D<double>>& interpolator) {
-    return new AutoDiffCostFunction<InterpolatedCostFunctor, 1, 1>(
-        new InterpolatedCostFunctor(interpolator));
+  static ceres::CostFunction* Create(const Interpolator& interpolator) {
+    return new ceres::AutoDiffCostFunction<InterpolatedCostFunctor, 1, 1>(
+        interpolator);
   }
 
  private:
-  const CubicInterpolator<Grid1D<double>>& interpolator_;
+  const Interpolator& interpolator;
 };
 
 int main(int argc, char** argv) {
@@ -76,18 +68,19 @@
     values[i] = (i - 4.5) * (i - 4.5);
   }
 
-  Grid1D<double> array(values, 0, kNumSamples);
-  CubicInterpolator<Grid1D<double>> interpolator(array);
+  ceres::Grid1D<double> array(values, 0, kNumSamples);
+  Interpolator interpolator(array);
 
   double x = 1.0;
-  Problem problem;
-  CostFunction* cost_function = InterpolatedCostFunctor::Create(interpolator);
-  problem.AddResidualBlock(cost_function, NULL, &x);
+  ceres::Problem problem;
+  ceres::CostFunction* cost_function =
+      InterpolatedCostFunctor::Create(interpolator);
+  problem.AddResidualBlock(cost_function, nullptr, &x);
 
-  Solver::Options options;
+  ceres::Solver::Options options;
   options.minimizer_progress_to_stdout = true;
-  Solver::Summary summary;
-  Solve(options, &problem, &summary);
+  ceres::Solver::Summary summary;
+  ceres::Solve(options, &problem, &summary);
   std::cout << summary.BriefReport() << "\n";
   std::cout << "Expected x: 4.5. Actual x : " << x << std::endl;
   return 0;
diff --git a/examples/simple_bundle_adjuster.cc b/examples/simple_bundle_adjuster.cc
index 8180d73..bb0ba1c 100644
--- a/examples/simple_bundle_adjuster.cc
+++ b/examples/simple_bundle_adjuster.cc
@@ -1,5 +1,5 @@
 // Ceres Solver - A fast non-linear least squares minimizer
-// Copyright 2015 Google Inc. All rights reserved.
+// Copyright 2023 Google Inc. All rights reserved.
 // http://ceres-solver.org/
 //
 // Redistribution and use in source and binary forms, with or without
@@ -66,7 +66,7 @@
 
   bool LoadFile(const char* filename) {
     FILE* fptr = fopen(filename, "r");
-    if (fptr == NULL) {
+    if (fptr == nullptr) {
       return false;
     };
 
@@ -164,8 +164,8 @@
   // the client code.
   static ceres::CostFunction* Create(const double observed_x,
                                      const double observed_y) {
-    return (new ceres::AutoDiffCostFunction<SnavelyReprojectionError, 2, 9, 3>(
-        new SnavelyReprojectionError(observed_x, observed_y)));
+    return new ceres::AutoDiffCostFunction<SnavelyReprojectionError, 2, 9, 3>(
+        observed_x, observed_y);
   }
 
   double observed_x;
@@ -198,7 +198,7 @@
     ceres::CostFunction* cost_function = SnavelyReprojectionError::Create(
         observations[2 * i + 0], observations[2 * i + 1]);
     problem.AddResidualBlock(cost_function,
-                             NULL /* squared loss */,
+                             nullptr /* squared loss */,
                              bal_problem.mutable_camera_for_observation(i),
                              bal_problem.mutable_point_for_observation(i));
   }
diff --git a/examples/slam/CMakeLists.txt b/examples/slam/CMakeLists.txt
index c72aa16..d0b03d7 100644
--- a/examples/slam/CMakeLists.txt
+++ b/examples/slam/CMakeLists.txt
@@ -1,5 +1,5 @@
 # Ceres Solver - A fast non-linear least squares minimizer
-# Copyright 2016 Google Inc. All rights reserved.
+# Copyright 2023 Google Inc. All rights reserved.
 # http://ceres-solver.org/
 #
 # Redistribution and use in source and binary forms, with or without
diff --git a/examples/slam/common/read_g2o.h b/examples/slam/common/read_g2o.h
index fea32e9..490b054 100644
--- a/examples/slam/common/read_g2o.h
+++ b/examples/slam/common/read_g2o.h
@@ -1,5 +1,5 @@
 // Ceres Solver - A fast non-linear least squares minimizer
-// Copyright 2016 Google Inc. All rights reserved.
+// Copyright 2023 Google Inc. All rights reserved.
 // http://ceres-solver.org/
 //
 // Redistribution and use in source and binary forms, with or without
@@ -38,8 +38,7 @@
 
 #include "glog/logging.h"
 
-namespace ceres {
-namespace examples {
+namespace ceres::examples {
 
 // Reads a single pose from the input and inserts it into the map. Returns false
 // if there is a duplicate entry.
@@ -60,7 +59,7 @@
   return true;
 }
 
-// Reads the contraints between two vertices in the pose graph
+// Reads the constraints between two vertices in the pose graph
 template <typename Constraint, typename Allocator>
 void ReadConstraint(std::ifstream* infile,
                     std::vector<Constraint, Allocator>* constraints) {
@@ -104,8 +103,8 @@
 bool ReadG2oFile(const std::string& filename,
                  std::map<int, Pose, std::less<int>, MapAllocator>* poses,
                  std::vector<Constraint, VectorAllocator>* constraints) {
-  CHECK(poses != NULL);
-  CHECK(constraints != NULL);
+  CHECK(poses != nullptr);
+  CHECK(constraints != nullptr);
 
   poses->clear();
   constraints->clear();
@@ -137,7 +136,6 @@
   return true;
 }
 
-}  // namespace examples
-}  // namespace ceres
+}  // namespace ceres::examples
 
 #endif  // EXAMPLES_CERES_READ_G2O_H_
diff --git a/examples/slam/pose_graph_2d/CMakeLists.txt b/examples/slam/pose_graph_2d/CMakeLists.txt
index 20af056..87943ec 100644
--- a/examples/slam/pose_graph_2d/CMakeLists.txt
+++ b/examples/slam/pose_graph_2d/CMakeLists.txt
@@ -1,5 +1,5 @@
 # Ceres Solver - A fast non-linear least squares minimizer
-# Copyright 2016 Google Inc. All rights reserved.
+# Copyright 2023 Google Inc. All rights reserved.
 # http://ceres-solver.org/
 #
 # Redistribution and use in source and binary forms, with or without
@@ -30,10 +30,10 @@
 
 if (GFLAGS)
   add_executable(pose_graph_2d
-    angle_local_parameterization.h
+    angle_manifold.h
     normalize_angle.h
     pose_graph_2d.cc
     pose_graph_2d_error_term.h
     types.h)
-  target_link_libraries(pose_graph_2d Ceres::ceres gflags)
+  target_link_libraries(pose_graph_2d PRIVATE Ceres::ceres gflags)
 endif (GFLAGS)
diff --git a/examples/slam/pose_graph_2d/angle_local_parameterization.h b/examples/slam/pose_graph_2d/angle_manifold.h
similarity index 63%
rename from examples/slam/pose_graph_2d/angle_local_parameterization.h
rename to examples/slam/pose_graph_2d/angle_manifold.h
index a81637c..456d923 100644
--- a/examples/slam/pose_graph_2d/angle_local_parameterization.h
+++ b/examples/slam/pose_graph_2d/angle_manifold.h
@@ -1,5 +1,5 @@
 // Ceres Solver - A fast non-linear least squares minimizer
-// Copyright 2016 Google Inc. All rights reserved.
+// Copyright 2023 Google Inc. All rights reserved.
 // http://ceres-solver.org/
 //
 // Redistribution and use in source and binary forms, with or without
@@ -27,38 +27,42 @@
 // POSSIBILITY OF SUCH DAMAGE.
 //
 // Author: vitus@google.com (Michael Vitus)
+//         sameeragarwal@google.com (Sameer Agarwal)
 
-#ifndef CERES_EXAMPLES_POSE_GRAPH_2D_ANGLE_LOCAL_PARAMETERIZATION_H_
-#define CERES_EXAMPLES_POSE_GRAPH_2D_ANGLE_LOCAL_PARAMETERIZATION_H_
+#ifndef CERES_EXAMPLES_POSE_GRAPH_2D_ANGLE_MANIFOLD_H_
+#define CERES_EXAMPLES_POSE_GRAPH_2D_ANGLE_MANIFOLD_H_
 
-#include "ceres/local_parameterization.h"
+#include "ceres/autodiff_manifold.h"
 #include "normalize_angle.h"
 
-namespace ceres {
-namespace examples {
+namespace ceres::examples {
 
-// Defines a local parameterization for updating the angle to be constrained in
-// [-pi to pi).
-class AngleLocalParameterization {
+// Defines a manifold for updating the angle to be constrained in [-pi to pi).
+class AngleManifold {
  public:
   template <typename T>
-  bool operator()(const T* theta_radians,
-                  const T* delta_theta_radians,
-                  T* theta_radians_plus_delta) const {
-    *theta_radians_plus_delta =
-        NormalizeAngle(*theta_radians + *delta_theta_radians);
+  bool Plus(const T* x_radians,
+            const T* delta_radians,
+            T* x_plus_delta_radians) const {
+    *x_plus_delta_radians = NormalizeAngle(*x_radians + *delta_radians);
+    return true;
+  }
+
+  template <typename T>
+  bool Minus(const T* y_radians,
+             const T* x_radians,
+             T* y_minus_x_radians) const {
+    *y_minus_x_radians =
+        NormalizeAngle(*y_radians) - NormalizeAngle(*x_radians);
 
     return true;
   }
 
-  static ceres::LocalParameterization* Create() {
-    return (new ceres::AutoDiffLocalParameterization<AngleLocalParameterization,
-                                                     1,
-                                                     1>);
+  static ceres::Manifold* Create() {
+    return new ceres::AutoDiffManifold<AngleManifold, 1, 1>;
   }
 };
 
-}  // namespace examples
-}  // namespace ceres
+}  // namespace ceres::examples
 
-#endif  // CERES_EXAMPLES_POSE_GRAPH_2D_ANGLE_LOCAL_PARAMETERIZATION_H_
+#endif  // CERES_EXAMPLES_POSE_GRAPH_2D_ANGLE_MANIFOLD_H_
diff --git a/examples/slam/pose_graph_2d/normalize_angle.h b/examples/slam/pose_graph_2d/normalize_angle.h
index c215671..0602878 100644
--- a/examples/slam/pose_graph_2d/normalize_angle.h
+++ b/examples/slam/pose_graph_2d/normalize_angle.h
@@ -1,5 +1,5 @@
 // Ceres Solver - A fast non-linear least squares minimizer
-// Copyright 2016 Google Inc. All rights reserved.
+// Copyright 2023 Google Inc. All rights reserved.
 // http://ceres-solver.org/
 //
 // Redistribution and use in source and binary forms, with or without
@@ -35,19 +35,17 @@
 
 #include "ceres/ceres.h"
 
-namespace ceres {
-namespace examples {
+namespace ceres::examples {
 
 // Normalizes the angle in radians between [-pi and pi).
 template <typename T>
 inline T NormalizeAngle(const T& angle_radians) {
   // Use ceres::floor because it is specialized for double and Jet types.
-  T two_pi(2.0 * M_PI);
+  T two_pi(2.0 * constants::pi);
   return angle_radians -
-         two_pi * ceres::floor((angle_radians + T(M_PI)) / two_pi);
+         two_pi * ceres::floor((angle_radians + T(constants::pi)) / two_pi);
 }
 
-}  // namespace examples
-}  // namespace ceres
+}  // namespace ceres::examples
 
 #endif  // CERES_EXAMPLES_POSE_GRAPH_2D_NORMALIZE_ANGLE_H_
diff --git a/examples/slam/pose_graph_2d/pose_graph_2d.cc b/examples/slam/pose_graph_2d/pose_graph_2d.cc
index 1172123..3ebae3f 100644
--- a/examples/slam/pose_graph_2d/pose_graph_2d.cc
+++ b/examples/slam/pose_graph_2d/pose_graph_2d.cc
@@ -1,5 +1,5 @@
 // Ceres Solver - A fast non-linear least squares minimizer
-// Copyright 2016 Google Inc. All rights reserved.
+// Copyright 2023 Google Inc. All rights reserved.
 // http://ceres-solver.org/
 //
 // Redistribution and use in source and binary forms, with or without
@@ -39,7 +39,7 @@
 #include <string>
 #include <vector>
 
-#include "angle_local_parameterization.h"
+#include "angle_manifold.h"
 #include "ceres/ceres.h"
 #include "common/read_g2o.h"
 #include "gflags/gflags.h"
@@ -49,8 +49,7 @@
 
 DEFINE_string(input, "", "The pose graph definition filename in g2o format.");
 
-namespace ceres {
-namespace examples {
+namespace ceres::examples {
 namespace {
 
 // Constructs the nonlinear least squares optimization problem from the pose
@@ -58,29 +57,21 @@
 void BuildOptimizationProblem(const std::vector<Constraint2d>& constraints,
                               std::map<int, Pose2d>* poses,
                               ceres::Problem* problem) {
-  CHECK(poses != NULL);
-  CHECK(problem != NULL);
+  CHECK(poses != nullptr);
+  CHECK(problem != nullptr);
   if (constraints.empty()) {
     LOG(INFO) << "No constraints, no problem to optimize.";
     return;
   }
 
-  ceres::LossFunction* loss_function = NULL;
-  ceres::LocalParameterization* angle_local_parameterization =
-      AngleLocalParameterization::Create();
+  ceres::LossFunction* loss_function = nullptr;
+  ceres::Manifold* angle_manifold = AngleManifold::Create();
 
-  for (std::vector<Constraint2d>::const_iterator constraints_iter =
-           constraints.begin();
-       constraints_iter != constraints.end();
-       ++constraints_iter) {
-    const Constraint2d& constraint = *constraints_iter;
-
-    std::map<int, Pose2d>::iterator pose_begin_iter =
-        poses->find(constraint.id_begin);
+  for (const auto& constraint : constraints) {
+    auto pose_begin_iter = poses->find(constraint.id_begin);
     CHECK(pose_begin_iter != poses->end())
         << "Pose with ID: " << constraint.id_begin << " not found.";
-    std::map<int, Pose2d>::iterator pose_end_iter =
-        poses->find(constraint.id_end);
+    auto pose_end_iter = poses->find(constraint.id_end);
     CHECK(pose_end_iter != poses->end())
         << "Pose with ID: " << constraint.id_end << " not found.";
 
@@ -98,10 +89,8 @@
                               &pose_end_iter->second.y,
                               &pose_end_iter->second.yaw_radians);
 
-    problem->SetParameterization(&pose_begin_iter->second.yaw_radians,
-                                 angle_local_parameterization);
-    problem->SetParameterization(&pose_end_iter->second.yaw_radians,
-                                 angle_local_parameterization);
+    problem->SetManifold(&pose_begin_iter->second.yaw_radians, angle_manifold);
+    problem->SetManifold(&pose_end_iter->second.yaw_radians, angle_manifold);
   }
 
   // The pose graph optimization problem has three DOFs that are not fully
@@ -111,7 +100,7 @@
   // internal damping which mitigate this issue, but it is better to properly
   // constrain the gauge freedom. This can be done by setting one of the poses
   // as constant so the optimizer cannot change it.
-  std::map<int, Pose2d>::iterator pose_start_iter = poses->begin();
+  auto pose_start_iter = poses->begin();
   CHECK(pose_start_iter != poses->end()) << "There are no poses.";
   problem->SetParameterBlockConstant(&pose_start_iter->second.x);
   problem->SetParameterBlockConstant(&pose_start_iter->second.y);
@@ -120,7 +109,7 @@
 
 // Returns true if the solve was successful.
 bool SolveOptimizationProblem(ceres::Problem* problem) {
-  CHECK(problem != NULL);
+  CHECK(problem != nullptr);
 
   ceres::Solver::Options options;
   options.max_num_iterations = 100;
@@ -143,10 +132,7 @@
     std::cerr << "Error opening the file: " << filename << '\n';
     return false;
   }
-  for (std::map<int, Pose2d>::const_iterator poses_iter = poses.begin();
-       poses_iter != poses.end();
-       ++poses_iter) {
-    const std::map<int, Pose2d>::value_type& pair = *poses_iter;
+  for (const auto& pair : poses) {
     outfile << pair.first << " " << pair.second.x << " " << pair.second.y << ' '
             << pair.second.yaw_radians << '\n';
   }
@@ -154,8 +140,7 @@
 }
 
 }  // namespace
-}  // namespace examples
-}  // namespace ceres
+}  // namespace ceres::examples
 
 int main(int argc, char** argv) {
   google::InitGoogleLogging(argv[0]);
diff --git a/examples/slam/pose_graph_2d/pose_graph_2d_error_term.h b/examples/slam/pose_graph_2d/pose_graph_2d_error_term.h
index 2df31f6..3d34f8d 100644
--- a/examples/slam/pose_graph_2d/pose_graph_2d_error_term.h
+++ b/examples/slam/pose_graph_2d/pose_graph_2d_error_term.h
@@ -1,5 +1,5 @@
 // Ceres Solver - A fast non-linear least squares minimizer
-// Copyright 2016 Google Inc. All rights reserved.
+// Copyright 2024 Google Inc. All rights reserved.
 // http://ceres-solver.org/
 //
 // Redistribution and use in source and binary forms, with or without
@@ -34,9 +34,9 @@
 #define CERES_EXAMPLES_POSE_GRAPH_2D_POSE_GRAPH_2D_ERROR_TERM_H_
 
 #include "Eigen/Core"
+#include "ceres/autodiff_cost_function.h"
 
-namespace ceres {
-namespace examples {
+namespace ceres::examples {
 
 template <typename T>
 Eigen::Matrix<T, 2, 2> RotationMatrix2D(T yaw_radians) {
@@ -96,10 +96,9 @@
                                      double y_ab,
                                      double yaw_ab_radians,
                                      const Eigen::Matrix3d& sqrt_information) {
-    return (new ceres::
-                AutoDiffCostFunction<PoseGraph2dErrorTerm, 3, 1, 1, 1, 1, 1, 1>(
-                    new PoseGraph2dErrorTerm(
-                        x_ab, y_ab, yaw_ab_radians, sqrt_information)));
+    return new ceres::
+        AutoDiffCostFunction<PoseGraph2dErrorTerm, 3, 1, 1, 1, 1, 1, 1>(
+            x_ab, y_ab, yaw_ab_radians, sqrt_information);
   }
 
   EIGEN_MAKE_ALIGNED_OPERATOR_NEW
@@ -113,7 +112,6 @@
   const Eigen::Matrix3d sqrt_information_;
 };
 
-}  // namespace examples
-}  // namespace ceres
+}  // namespace ceres::examples
 
 #endif  // CERES_EXAMPLES_POSE_GRAPH_2D_POSE_GRAPH_2D_ERROR_TERM_H_
diff --git a/examples/slam/pose_graph_2d/types.h b/examples/slam/pose_graph_2d/types.h
index 3c13824..caf4ccb 100644
--- a/examples/slam/pose_graph_2d/types.h
+++ b/examples/slam/pose_graph_2d/types.h
@@ -1,5 +1,5 @@
 // Ceres Solver - A fast non-linear least squares minimizer
-// Copyright 2016 Google Inc. All rights reserved.
+// Copyright 2023 Google Inc. All rights reserved.
 // http://ceres-solver.org/
 //
 // Redistribution and use in source and binary forms, with or without
@@ -40,8 +40,7 @@
 #include "Eigen/Core"
 #include "normalize_angle.h"
 
-namespace ceres {
-namespace examples {
+namespace ceres::examples {
 
 // The state for each vertex in the pose graph.
 struct Pose2d {
@@ -95,7 +94,6 @@
   return input;
 }
 
-}  // namespace examples
-}  // namespace ceres
+}  // namespace ceres::examples
 
 #endif  // CERES_EXAMPLES_POSE_GRAPH_2D_TYPES_H_
diff --git a/examples/slam/pose_graph_3d/CMakeLists.txt b/examples/slam/pose_graph_3d/CMakeLists.txt
index b6421cc..544b00e 100644
--- a/examples/slam/pose_graph_3d/CMakeLists.txt
+++ b/examples/slam/pose_graph_3d/CMakeLists.txt
@@ -1,5 +1,5 @@
 # Ceres Solver - A fast non-linear least squares minimizer
-# Copyright 2016 Google Inc. All rights reserved.
+# Copyright 2023 Google Inc. All rights reserved.
 # http://ceres-solver.org/
 #
 # Redistribution and use in source and binary forms, with or without
@@ -30,5 +30,5 @@
 
 if (GFLAGS)
   add_executable(pose_graph_3d pose_graph_3d.cc)
-  target_link_libraries(pose_graph_3d Ceres::ceres gflags)
+  target_link_libraries(pose_graph_3d PRIVATE Ceres::ceres gflags)
 endif (GFLAGS)
diff --git a/examples/slam/pose_graph_3d/pose_graph_3d.cc b/examples/slam/pose_graph_3d/pose_graph_3d.cc
index 2f8d6a4..522e2a1 100644
--- a/examples/slam/pose_graph_3d/pose_graph_3d.cc
+++ b/examples/slam/pose_graph_3d/pose_graph_3d.cc
@@ -1,5 +1,5 @@
 // Ceres Solver - A fast non-linear least squares minimizer
-// Copyright 2016 Google Inc. All rights reserved.
+// Copyright 2023 Google Inc. All rights reserved.
 // http://ceres-solver.org/
 //
 // Redistribution and use in source and binary forms, with or without
@@ -41,8 +41,7 @@
 
 DEFINE_string(input, "", "The pose graph definition filename in g2o format.");
 
-namespace ceres {
-namespace examples {
+namespace ceres::examples {
 namespace {
 
 // Constructs the nonlinear least squares optimization problem from the pose
@@ -50,27 +49,21 @@
 void BuildOptimizationProblem(const VectorOfConstraints& constraints,
                               MapOfPoses* poses,
                               ceres::Problem* problem) {
-  CHECK(poses != NULL);
-  CHECK(problem != NULL);
+  CHECK(poses != nullptr);
+  CHECK(problem != nullptr);
   if (constraints.empty()) {
     LOG(INFO) << "No constraints, no problem to optimize.";
     return;
   }
 
-  ceres::LossFunction* loss_function = NULL;
-  ceres::LocalParameterization* quaternion_local_parameterization =
-      new EigenQuaternionParameterization;
+  ceres::LossFunction* loss_function = nullptr;
+  ceres::Manifold* quaternion_manifold = new EigenQuaternionManifold;
 
-  for (VectorOfConstraints::const_iterator constraints_iter =
-           constraints.begin();
-       constraints_iter != constraints.end();
-       ++constraints_iter) {
-    const Constraint3d& constraint = *constraints_iter;
-
-    MapOfPoses::iterator pose_begin_iter = poses->find(constraint.id_begin);
+  for (const auto& constraint : constraints) {
+    auto pose_begin_iter = poses->find(constraint.id_begin);
     CHECK(pose_begin_iter != poses->end())
         << "Pose with ID: " << constraint.id_begin << " not found.";
-    MapOfPoses::iterator pose_end_iter = poses->find(constraint.id_end);
+    auto pose_end_iter = poses->find(constraint.id_end);
     CHECK(pose_end_iter != poses->end())
         << "Pose with ID: " << constraint.id_end << " not found.";
 
@@ -87,10 +80,10 @@
                               pose_end_iter->second.p.data(),
                               pose_end_iter->second.q.coeffs().data());
 
-    problem->SetParameterization(pose_begin_iter->second.q.coeffs().data(),
-                                 quaternion_local_parameterization);
-    problem->SetParameterization(pose_end_iter->second.q.coeffs().data(),
-                                 quaternion_local_parameterization);
+    problem->SetManifold(pose_begin_iter->second.q.coeffs().data(),
+                         quaternion_manifold);
+    problem->SetManifold(pose_end_iter->second.q.coeffs().data(),
+                         quaternion_manifold);
   }
 
   // The pose graph optimization problem has six DOFs that are not fully
@@ -100,7 +93,7 @@
   // internal damping which mitigates this issue, but it is better to properly
   // constrain the gauge freedom. This can be done by setting one of the poses
   // as constant so the optimizer cannot change it.
-  MapOfPoses::iterator pose_start_iter = poses->begin();
+  auto pose_start_iter = poses->begin();
   CHECK(pose_start_iter != poses->end()) << "There are no poses.";
   problem->SetParameterBlockConstant(pose_start_iter->second.p.data());
   problem->SetParameterBlockConstant(pose_start_iter->second.q.coeffs().data());
@@ -108,7 +101,7 @@
 
 // Returns true if the solve was successful.
 bool SolveOptimizationProblem(ceres::Problem* problem) {
-  CHECK(problem != NULL);
+  CHECK(problem != nullptr);
 
   ceres::Solver::Options options;
   options.max_num_iterations = 200;
@@ -130,18 +123,7 @@
     LOG(ERROR) << "Error opening the file: " << filename;
     return false;
   }
-  for (std::map<int,
-                Pose3d,
-                std::less<int>,
-                Eigen::aligned_allocator<std::pair<const int, Pose3d>>>::
-           const_iterator poses_iter = poses.begin();
-       poses_iter != poses.end();
-       ++poses_iter) {
-    const std::map<int,
-                   Pose3d,
-                   std::less<int>,
-                   Eigen::aligned_allocator<std::pair<const int, Pose3d>>>::
-        value_type& pair = *poses_iter;
+  for (const auto& pair : poses) {
     outfile << pair.first << " " << pair.second.p.transpose() << " "
             << pair.second.q.x() << " " << pair.second.q.y() << " "
             << pair.second.q.z() << " " << pair.second.q.w() << '\n';
@@ -150,8 +132,7 @@
 }
 
 }  // namespace
-}  // namespace examples
-}  // namespace ceres
+}  // namespace ceres::examples
 
 int main(int argc, char** argv) {
   google::InitGoogleLogging(argv[0]);
diff --git a/examples/slam/pose_graph_3d/pose_graph_3d_error_term.h b/examples/slam/pose_graph_3d/pose_graph_3d_error_term.h
index 1f3e8de..b1c0138 100644
--- a/examples/slam/pose_graph_3d/pose_graph_3d_error_term.h
+++ b/examples/slam/pose_graph_3d/pose_graph_3d_error_term.h
@@ -1,5 +1,5 @@
 // Ceres Solver - A fast non-linear least squares minimizer
-// Copyright 2016 Google Inc. All rights reserved.
+// Copyright 2023 Google Inc. All rights reserved.
 // http://ceres-solver.org/
 //
 // Redistribution and use in source and binary forms, with or without
@@ -31,12 +31,13 @@
 #ifndef EXAMPLES_CERES_POSE_GRAPH_3D_ERROR_TERM_H_
 #define EXAMPLES_CERES_POSE_GRAPH_3D_ERROR_TERM_H_
 
+#include <utility>
+
 #include "Eigen/Core"
 #include "ceres/autodiff_cost_function.h"
 #include "types.h"
 
-namespace ceres {
-namespace examples {
+namespace ceres::examples {
 
 // Computes the error term for two poses that have a relative pose measurement
 // between them. Let the hat variables be the measurement. We have two poses x_a
@@ -69,9 +70,10 @@
 // where I is the information matrix which is the inverse of the covariance.
 class PoseGraph3dErrorTerm {
  public:
-  PoseGraph3dErrorTerm(const Pose3d& t_ab_measured,
-                       const Eigen::Matrix<double, 6, 6>& sqrt_information)
-      : t_ab_measured_(t_ab_measured), sqrt_information_(sqrt_information) {}
+  PoseGraph3dErrorTerm(Pose3d t_ab_measured,
+                       Eigen::Matrix<double, 6, 6> sqrt_information)
+      : t_ab_measured_(std::move(t_ab_measured)),
+        sqrt_information_(std::move(sqrt_information)) {}
 
   template <typename T>
   bool operator()(const T* const p_a_ptr,
@@ -114,7 +116,7 @@
       const Pose3d& t_ab_measured,
       const Eigen::Matrix<double, 6, 6>& sqrt_information) {
     return new ceres::AutoDiffCostFunction<PoseGraph3dErrorTerm, 6, 3, 4, 3, 4>(
-        new PoseGraph3dErrorTerm(t_ab_measured, sqrt_information));
+        t_ab_measured, sqrt_information);
   }
 
   EIGEN_MAKE_ALIGNED_OPERATOR_NEW
@@ -126,7 +128,6 @@
   const Eigen::Matrix<double, 6, 6> sqrt_information_;
 };
 
-}  // namespace examples
-}  // namespace ceres
+}  // namespace ceres::examples
 
 #endif  // EXAMPLES_CERES_POSE_GRAPH_3D_ERROR_TERM_H_
diff --git a/examples/slam/pose_graph_3d/types.h b/examples/slam/pose_graph_3d/types.h
index d3f19ed..207dd5d 100644
--- a/examples/slam/pose_graph_3d/types.h
+++ b/examples/slam/pose_graph_3d/types.h
@@ -1,5 +1,5 @@
 // Ceres Solver - A fast non-linear least squares minimizer
-// Copyright 2016 Google Inc. All rights reserved.
+// Copyright 2023 Google Inc. All rights reserved.
 // http://ceres-solver.org/
 //
 // Redistribution and use in source and binary forms, with or without
@@ -31,6 +31,7 @@
 #ifndef EXAMPLES_CERES_TYPES_H_
 #define EXAMPLES_CERES_TYPES_H_
 
+#include <functional>
 #include <istream>
 #include <map>
 #include <string>
@@ -39,8 +40,7 @@
 #include "Eigen/Core"
 #include "Eigen/Geometry"
 
-namespace ceres {
-namespace examples {
+namespace ceres::examples {
 
 struct Pose3d {
   Eigen::Vector3d p;
@@ -61,11 +61,11 @@
   return input;
 }
 
-typedef std::map<int,
-                 Pose3d,
-                 std::less<int>,
-                 Eigen::aligned_allocator<std::pair<const int, Pose3d>>>
-    MapOfPoses;
+using MapOfPoses =
+    std::map<int,
+             Pose3d,
+             std::less<int>,
+             Eigen::aligned_allocator<std::pair<const int, Pose3d>>>;
 
 // The constraint between two vertices in the pose graph. The constraint is the
 // transformation from vertex id_begin to vertex id_end.
@@ -103,10 +103,9 @@
   return input;
 }
 
-typedef std::vector<Constraint3d, Eigen::aligned_allocator<Constraint3d>>
-    VectorOfConstraints;
+using VectorOfConstraints =
+    std::vector<Constraint3d, Eigen::aligned_allocator<Constraint3d>>;
 
-}  // namespace examples
-}  // namespace ceres
+}  // namespace ceres::examples
 
 #endif  // EXAMPLES_CERES_TYPES_H_
diff --git a/examples/snavely_reprojection_error.h b/examples/snavely_reprojection_error.h
index eb39d23..aaf0c6c 100644
--- a/examples/snavely_reprojection_error.h
+++ b/examples/snavely_reprojection_error.h
@@ -1,5 +1,5 @@
 // Ceres Solver - A fast non-linear least squares minimizer
-// Copyright 2015 Google Inc. All rights reserved.
+// Copyright 2024 Google Inc. All rights reserved.
 // http://ceres-solver.org/
 //
 // Redistribution and use in source and binary forms, with or without
@@ -41,10 +41,10 @@
 #ifndef CERES_EXAMPLES_SNAVELY_REPROJECTION_ERROR_H_
 #define CERES_EXAMPLES_SNAVELY_REPROJECTION_ERROR_H_
 
+#include "ceres/autodiff_cost_function.h"
 #include "ceres/rotation.h"
 
-namespace ceres {
-namespace examples {
+namespace ceres::examples {
 
 // Templated pinhole camera model for used with Ceres.  The camera is
 // parameterized using 9 parameters: 3 for rotation, 3 for translation, 1 for
@@ -95,8 +95,8 @@
   // the client code.
   static ceres::CostFunction* Create(const double observed_x,
                                      const double observed_y) {
-    return (new ceres::AutoDiffCostFunction<SnavelyReprojectionError, 2, 9, 3>(
-        new SnavelyReprojectionError(observed_x, observed_y)));
+    return new ceres::AutoDiffCostFunction<SnavelyReprojectionError, 2, 9, 3>(
+        observed_x, observed_y);
   }
 
   double observed_x;
@@ -123,7 +123,7 @@
     // We use QuaternionRotatePoint as it does not assume that the
     // quaternion is normalized, since one of the ways to run the
     // bundle adjuster is to let Ceres optimize all 4 quaternion
-    // parameters without a local parameterization.
+    // parameters without using a Quaternion manifold.
     T p[3];
     QuaternionRotatePoint(camera, point, p);
 
@@ -160,20 +160,15 @@
   // the client code.
   static ceres::CostFunction* Create(const double observed_x,
                                      const double observed_y) {
-    return (
-        new ceres::AutoDiffCostFunction<SnavelyReprojectionErrorWithQuaternions,
-                                        2,
-                                        10,
-                                        3>(
-            new SnavelyReprojectionErrorWithQuaternions(observed_x,
-                                                        observed_y)));
+    return new ceres::
+        AutoDiffCostFunction<SnavelyReprojectionErrorWithQuaternions, 2, 10, 3>(
+            observed_x, observed_y);
   }
 
   double observed_x;
   double observed_y;
 };
 
-}  // namespace examples
-}  // namespace ceres
+}  // namespace ceres::examples
 
 #endif  // CERES_EXAMPLES_SNAVELY_REPROJECTION_ERROR_H_