Struct ParallelCsr32Operator 

pub struct ParallelCsr32Operator<'a> { /* private fields */ }

Read-only parallel CSR operator wrapper. It shares the validated CSR storage and only changes the row execution policy used by apply.

Implementations

impl<'a> ParallelCsr32Operator<'a>

pub fn new(matrix: &'a Csr32Matrix) -> ParallelCsr32Operator<'a>

Examples found in repository

examples/fem_structural_auto.rs (line 389)

fn main() -> Result<(), Box<dyn Error>> {
    let args = Args::parse()?;
    println!(
        "HyBIT {} structural-auto FEM benchmark",
        env!("CARGO_PKG_VERSION")
    );
    println!("matrix             : {}", args.matrix.display());
    println!("coordinates        : {}", args.coordinates.display());

    let load_start = Instant::now();
    let (matrix, mm) = read_matrix_market(&args.matrix)?;
    let matrix_load_seconds = load_start.elapsed().as_secs_f64();
    let coord_start = Instant::now();
    let coordinates = read_coordinates(&args.coordinates)?;
    let coord_load_seconds = coord_start.elapsed().as_secs_f64();
    let profile = analyze_csr32(&matrix)?;

    println!(
        "Matrix Market      : {:?}, {} input entries -> {} CSR nnz",
        mm.symmetry, mm.input_entries, mm.csr_nnz
    );
    println!("dimensions         : {} x {}", profile.nrows, profile.ncols);
    println!("nnz                : {}", profile.nnz);
    println!(
        "CSR storage        : {:.3} MiB",
        mib(matrix.storage_bytes())
    );
    println!("matrix load        : {:.3} ms", matrix_load_seconds * 1.0e3);
    println!("coordinate nodes   : {}", coordinates.len());
    println!("coordinate load    : {:.3} ms", coord_load_seconds * 1.0e3);

    if !profile.square || !profile.full_diagonal || !profile.positive_diagonal {
        return Err(
            "rigid-body two-level PCG requires a square matrix with a complete positive diagonal"
                .into(),
        );
    }
    if matrix.nrows() != coordinates.len() * 3 {
        return Err(format!(
            "matrix/coordinate mismatch: {} matrix rows != {} coordinate nodes * 3",
            matrix.nrows(),
            coordinates.len()
        )
        .into());
    }

    let generated_rhs = args.rhs.is_none();
    let b = if let Some(path) = args.rhs.as_deref() {
        println!("RHS                : {}", path.display());
        load_rhs(path, matrix.nrows())?
    } else {
        println!("RHS                : generated as b=A*1 (known exact solution)");
        let ones = vec![1.0; matrix.ncols()];
        matrix.spmv(&ones)?
    };
    let mut solver = HybitSolver::new();
    solver.set_options(SolverOptions {
        relative_tolerance: args.relative_tolerance,
        absolute_tolerance: 0.0,
        max_iterations: args.max_iterations,
    })?;
    solver.set_structural_options(StructuralOptions {
        target_coarse_dimension: args.target_coarse_dimension,
        aggregation: args.aggregation,
        spmv_policy: args.spmv_policy,
        preconditioner_policy: args.preconditioner_policy,
        pcg_vector_policy: args.pcg_vector_policy,
    })?;

    let analysis = solver.analyze_csr32(&matrix)?;
    let mut prepared = solver.prepare_structural_csr32(&matrix, &analysis, &coordinates)?;
    println!("policy             : StructuralAuto/RigidBodyTwoLevel");
    println!("aggregation        : {:?}", prepared.aggregation());
    println!("SpMV policy        : {:?}", prepared.spmv_policy());
    println!(
        "precond policy     : {:?}",
        prepared.structural_preconditioner_policy()
    );
    println!("PCG vector policy  : {:?}", prepared.pcg_vector_policy());
    if prepared.parallel_spmv_enabled()
        || prepared.parallel_preconditioner_enabled()
        || prepared.parallel_pcg_vectors_enabled()
    {
        println!(
            "Rayon threads      : {}",
            ParallelCsr32Operator::new(&matrix).rayon_threads()
        );
    }
    if prepared.parallel_preconditioner_enabled() {
        println!(
            "parallel index     : {:.3} MiB",
            mib(prepared.parallel_preconditioner_index_bytes())
        );
    }
    println!("target coarse dim  : {}", args.target_coarse_dimension);
    println!(
        "aggregate nodes    : {} (auto-selected)",
        prepared.aggregate_nodes()
    );
    println!("fine block size    : 3");
    println!("aggregate count    : {}", prepared.aggregate_count());
    println!(
        "aggregate min/max  : {} / {} nodes",
        prepared.min_aggregate_nodes(),
        prepared.max_aggregate_nodes()
    );
    println!("modes/aggregate    : 6");
    println!("coarse dimension   : {}", prepared.coarse_dimension());
    println!(
        "base factor        : {:.3} MiB",
        mib(prepared.base_factor_bytes())
    );
    println!(
        "coarse factor      : {:.3} MiB",
        mib(prepared.coarse_factor_bytes())
    );
    println!(
        "geometry storage   : {:.3} MiB",
        mib(prepared.geometry_bytes())
    );
    println!(
        "total prec storage : {:.3} MiB",
        mib(prepared.preconditioner_bytes())
    );
    println!(
        "analysis           : {:.3} ms",
        prepared.analysis_seconds() * 1.0e3
    );
    println!(
        "prepare            : {:.3} ms",
        prepared.prepare_seconds() * 1.0e3
    );

    let mut x = vec![0.0; matrix.ncols()];
    let report = prepared.solve(&matrix, &b, &mut x)?;
    let verified = verified_relative_residual(&matrix, &b, &x)?;

    println!();
    println!("Structural Auto API: 3x3 Block-Jacobi + six rigid-body coarse modes");
    println!("status             : {:?}", report.status);
    println!("preconditioner     : {:?}", report.preconditioner);
    println!("iterations         : {}", report.iterations);
    println!("reported residual  : {:.6e}", report.relative_residual);
    println!("verified residual  : {:.6e}", verified);
    if generated_rhs {
        println!("relative x error   : {:.6e}", relative_error_to_ones(&x));
    }
    println!(
        "solve time         : {:.3} ms",
        report.solve_seconds * 1.0e3
    );
    println!(
        "total setup+solve  : {:.3} ms",
        (report.setup_seconds + report.solve_seconds) * 1.0e3
    );

    if !verified.is_finite() {
        return Err("non-finite independently verified residual".into());
    }
    Ok(())
}

examples/fem_structural_pcg_parallel.rs (line 324)

fn main() -> Result<(), Box<dyn Error>> {
    let args = Args::parse()?;
    println!(
        "HyBIT {} serial vs parallel/fused PCG vector benchmark",
        env!("CARGO_PKG_VERSION")
    );
    println!("matrix             : {}", args.matrix.display());
    println!("coordinates        : {}", args.coordinates.display());

    let load_start = Instant::now();
    let (matrix, mm) = read_matrix_market(&args.matrix)?;
    let matrix_load = load_start.elapsed().as_secs_f64();
    let coord_start = Instant::now();
    let coordinates = read_coordinates(&args.coordinates)?;
    let coord_load = coord_start.elapsed().as_secs_f64();
    let matrix_profile = analyze_csr32(&matrix)?;

    println!(
        "Matrix Market      : {:?}, {} input entries -> {} CSR nnz",
        mm.symmetry, mm.input_entries, mm.csr_nnz
    );
    println!(
        "dimensions         : {} x {}",
        matrix_profile.nrows, matrix_profile.ncols
    );
    println!("nnz                : {}", matrix_profile.nnz);
    println!(
        "CSR storage        : {:.3} MiB",
        mib(matrix.storage_bytes())
    );
    println!("matrix load        : {:.3} ms", matrix_load * 1.0e3);
    println!("coordinate nodes   : {}", coordinates.len());
    println!("coordinate load    : {:.3} ms", coord_load * 1.0e3);
    println!("aggregation        : {:?}", args.aggregation);
    println!("target coarse dim  : {}", args.target_coarse_dimension);

    if !matrix_profile.square || !matrix_profile.full_diagonal || !matrix_profile.positive_diagonal
    {
        return Err(
            "structural PCG benchmark requires a square matrix with a complete positive diagonal"
                .into(),
        );
    }
    if matrix.nrows() != coordinates.len() * 3 {
        return Err(format!(
            "matrix/coordinate mismatch: {} matrix rows != {} coordinate nodes * 3",
            matrix.nrows(),
            coordinates.len()
        )
        .into());
    }

    let b = if let Some(path) = args.rhs.as_deref() {
        println!("RHS                : {}", path.display());
        load_rhs(path, matrix.nrows())?
    } else {
        println!("RHS                : generated as b=A*1 (known exact solution)");
        matrix.spmv(&vec![1.0; matrix.ncols()])?
    };

    let aggregate_nodes =
        recommend_rigid_body_aggregate_nodes(coordinates.len(), args.target_coarse_dimension)?;
    let setup_start = Instant::now();
    let preconditioner = match args.aggregation {
        RigidBodyAggregation::Graph => {
            RigidBodyTwoLevelBlockJacobiPreconditioner::from_csr32_graph(
                &matrix,
                &coordinates,
                aggregate_nodes,
            )?
        }
        RigidBodyAggregation::Contiguous => RigidBodyTwoLevelBlockJacobiPreconditioner::from_csr32(
            &matrix,
            &coordinates,
            aggregate_nodes,
        )?,
        RigidBodyAggregation::Auto => unreachable!(),
    };
    let parallel_preconditioner = ParallelRigidBodyTwoLevelPreconditioner::new(&preconditioner)?;
    let parallel_operator = ParallelCsr32Operator::new(&matrix);
    let setup_seconds = setup_start.elapsed().as_secs_f64();

    println!("aggregate target   : {} nodes", aggregate_nodes);
    println!("aggregate count    : {}", preconditioner.aggregate_count());
    println!("coarse dimension   : {}", preconditioner.coarse_dimension());
    println!(
        "prec storage       : {:.3} MiB",
        mib(preconditioner.factor_bytes())
    );
    println!(
        "parallel index     : {:.3} MiB",
        mib(parallel_preconditioner.index_storage_bytes())
    );
    println!("setup              : {:.3} ms", setup_seconds * 1.0e3);
    println!(
        "Rayon threads      : {}",
        parallel_preconditioner.rayon_threads()
    );
    println!("vector chunk       : {} values", PARALLEL_PCG_VECTOR_CHUNK);

    let options = SolverOptions {
        relative_tolerance: args.relative_tolerance,
        absolute_tolerance: 0.0,
        max_iterations: args.max_iterations,
    };
    let b_norm = norm2(&b).max(f64::MIN_POSITIVE);

    let mut x_serial = vec![0.0; matrix.ncols()];
    let mut ws_serial = PcgWorkspace::new(matrix.nrows());
    let serial_start = Instant::now();
    let serial_out = pcg_with_workspace(
        &parallel_operator,
        &parallel_preconditioner,
        &b,
        &mut x_serial,
        options,
        &mut ws_serial,
    )?;
    let serial_wall = serial_start.elapsed().as_secs_f64();
    let serial_verified = verified_relative_residual(&matrix, &b, &x_serial)?;

    let mut x_parallel = vec![0.0; matrix.ncols()];
    let mut ws_parallel = PcgWorkspace::new(matrix.nrows());
    let parallel_start = Instant::now();
    let parallel_out = pcg_with_workspace_parallel_vectors(
        &parallel_operator,
        &parallel_preconditioner,
        &b,
        &mut x_parallel,
        options,
        &mut ws_parallel,
    )?;
    let parallel_wall = parallel_start.elapsed().as_secs_f64();
    let parallel_verified = verified_relative_residual(&matrix, &b, &x_parallel)?;

    println!();
    println!("[1/2] Production PCG vector kernels");
    println!("status             : {:?}", serial_out.status);
    println!("iterations         : {}", serial_out.iterations);
    println!(
        "reported residual  : {:.6e}",
        serial_out.final_residual / b_norm
    );
    println!("verified residual  : {:.6e}", serial_verified);
    println!("solve time         : {:.3} ms", serial_wall * 1.0e3);

    println!();
    println!("[2/2] Parallel/fused PCG vector kernels");
    println!("status             : {:?}", parallel_out.status);
    println!("iterations         : {}", parallel_out.iterations);
    println!(
        "reported residual  : {:.6e}",
        parallel_out.final_residual / b_norm
    );
    println!("verified residual  : {:.6e}", parallel_verified);
    println!("solve time         : {:.3} ms", parallel_wall * 1.0e3);

    println!();
    println!("Comparison");
    println!(
        "iteration ratio    : {:.3} (parallel vectors / production)",
        parallel_out.iterations as f64 / (serial_out.iterations.max(1) as f64)
    );
    println!(
        "solve-time ratio   : {:.3} (parallel vectors / production)",
        parallel_wall / serial_wall.max(f64::MIN_POSITIVE)
    );
    println!(
        "PCG speedup        : {:.3}x",
        serial_wall / parallel_wall.max(f64::MIN_POSITIVE)
    );
    println!(
        "solution delta     : {:.6e} relative L2",
        relative_difference(&x_serial, &x_parallel)
    );
    println!("note               : sparse operator and preconditioner are identical in both runs");

    Ok(())
}

examples/fem_structural_spmv.rs (line 319)

fn main() -> Result<(), Box<dyn Error>> {
    let args = Args::parse()?;
    println!(
        "HyBIT {} serial vs parallel CSR structural benchmark",
        env!("CARGO_PKG_VERSION")
    );
    println!("matrix             : {}", args.matrix.display());
    println!("coordinates        : {}", args.coordinates.display());

    let load_start = Instant::now();
    let (matrix, mm) = read_matrix_market(&args.matrix)?;
    let matrix_load = load_start.elapsed().as_secs_f64();
    let coord_start = Instant::now();
    let coordinates = read_coordinates(&args.coordinates)?;
    let coord_load = coord_start.elapsed().as_secs_f64();
    let profile = analyze_csr32(&matrix)?;

    println!(
        "Matrix Market      : {:?}, {} input entries -> {} CSR nnz",
        mm.symmetry, mm.input_entries, mm.csr_nnz
    );
    println!("dimensions         : {} x {}", profile.nrows, profile.ncols);
    println!("nnz                : {}", profile.nnz);
    println!(
        "CSR storage        : {:.3} MiB",
        mib(matrix.storage_bytes())
    );
    println!("matrix load        : {:.3} ms", matrix_load * 1.0e3);
    println!("coordinate nodes   : {}", coordinates.len());
    println!("coordinate load    : {:.3} ms", coord_load * 1.0e3);
    println!("aggregation        : {:?}", args.aggregation);
    println!("target coarse dim  : {}", args.target_coarse_dimension);
    println!("kernel repeats     : {}", args.kernel_repeats);

    if !profile.square || !profile.full_diagonal || !profile.positive_diagonal {
        return Err(
            "structural PCG benchmark requires a square matrix with a complete positive diagonal"
                .into(),
        );
    }
    if matrix.nrows() != coordinates.len() * 3 {
        return Err(format!(
            "matrix/coordinate mismatch: {} matrix rows != {} coordinate nodes * 3",
            matrix.nrows(),
            coordinates.len()
        )
        .into());
    }

    let b = if let Some(path) = args.rhs.as_deref() {
        println!("RHS                : {}", path.display());
        load_rhs(path, matrix.nrows())?
    } else {
        println!("RHS                : generated as b=A*1 (known exact solution)");
        matrix.spmv(&vec![1.0; matrix.ncols()])?
    };

    let aggregate_nodes =
        recommend_rigid_body_aggregate_nodes(coordinates.len(), args.target_coarse_dimension)?;
    let setup_start = Instant::now();
    let preconditioner = match args.aggregation {
        RigidBodyAggregation::Graph => {
            RigidBodyTwoLevelBlockJacobiPreconditioner::from_csr32_graph(
                &matrix,
                &coordinates,
                aggregate_nodes,
            )?
        }
        RigidBodyAggregation::Contiguous => RigidBodyTwoLevelBlockJacobiPreconditioner::from_csr32(
            &matrix,
            &coordinates,
            aggregate_nodes,
        )?,
        RigidBodyAggregation::Auto => unreachable!(),
    };
    let setup_seconds = setup_start.elapsed().as_secs_f64();
    let parallel = ParallelCsr32Operator::new(&matrix);

    println!("aggregate target   : {} nodes", aggregate_nodes);
    println!("aggregate count    : {}", preconditioner.aggregate_count());
    println!("coarse dimension   : {}", preconditioner.coarse_dimension());
    println!(
        "prec storage       : {:.3} MiB",
        mib(preconditioner.factor_bytes())
    );
    println!("setup              : {:.3} ms", setup_seconds * 1.0e3);
    println!("Rayon threads      : {}", parallel.rayon_threads());

    let mut ys = vec![0.0; matrix.nrows()];
    let mut yp = vec![0.0; matrix.nrows()];
    matrix.apply(&b, &mut ys)?;
    parallel.apply(&b, &mut yp)?;
    if ys != yp {
        return Err("serial and parallel CSR kernels produced different results".into());
    }

    let serial_start = Instant::now();
    for _ in 0..args.kernel_repeats {
        matrix.apply(&b, &mut ys)?;
    }
    let serial_spmv = serial_start.elapsed().as_secs_f64();

    let parallel_start = Instant::now();
    for _ in 0..args.kernel_repeats {
        parallel.apply(&b, &mut yp)?;
    }
    let parallel_spmv = parallel_start.elapsed().as_secs_f64();

    println!();
    println!("CSR SpMV microbenchmark");
    println!(
        "serial / call      : {:.3} ms",
        ms_per(serial_spmv, args.kernel_repeats)
    );
    println!(
        "parallel / call    : {:.3} ms",
        ms_per(parallel_spmv, args.kernel_repeats)
    );
    println!(
        "SpMV speedup       : {:.3}x",
        serial_spmv / parallel_spmv.max(f64::MIN_POSITIVE)
    );

    let options = SolverOptions {
        relative_tolerance: args.relative_tolerance,
        absolute_tolerance: 0.0,
        max_iterations: args.max_iterations,
    };

    let mut xs = vec![0.0; matrix.ncols()];
    let mut ws = PcgWorkspace::new(matrix.nrows());
    let serial_solve_start = Instant::now();
    let serial_out = pcg_with_workspace(&matrix, &preconditioner, &b, &mut xs, options, &mut ws)?;
    let serial_solve = serial_solve_start.elapsed().as_secs_f64();
    let serial_verified = verified_relative_residual(&matrix, &b, &xs)?;

    let mut xp = vec![0.0; matrix.ncols()];
    let mut wp = PcgWorkspace::new(matrix.nrows());
    let parallel_solve_start = Instant::now();
    let parallel_out =
        pcg_with_workspace(&parallel, &preconditioner, &b, &mut xp, options, &mut wp)?;
    let parallel_solve = parallel_solve_start.elapsed().as_secs_f64();
    let parallel_verified = verified_relative_residual(&matrix, &b, &xp)?;

    println!();
    println!("[1/2] Serial CSR PCG");
    println!("status             : {:?}", serial_out.status);
    println!("iterations         : {}", serial_out.iterations);
    println!(
        "reported residual  : {:.6e}",
        serial_out.final_residual / norm2(&b).max(f64::MIN_POSITIVE)
    );
    println!("verified residual  : {:.6e}", serial_verified);
    println!("solve time         : {:.3} ms", serial_solve * 1.0e3);

    println!();
    println!("[2/2] Parallel CSR PCG");
    println!("status             : {:?}", parallel_out.status);
    println!("iterations         : {}", parallel_out.iterations);
    println!(
        "reported residual  : {:.6e}",
        parallel_out.final_residual / norm2(&b).max(f64::MIN_POSITIVE)
    );
    println!("verified residual  : {:.6e}", parallel_verified);
    println!("solve time         : {:.3} ms", parallel_solve * 1.0e3);

    println!();
    println!("Comparison");
    println!(
        "iteration ratio    : {:.3} (parallel / serial)",
        parallel_out.iterations as f64 / serial_out.iterations.max(1) as f64
    );
    println!(
        "solve-time ratio   : {:.3} (parallel / serial)",
        parallel_solve / serial_solve.max(f64::MIN_POSITIVE)
    );
    println!(
        "PCG speedup        : {:.3}x",
        serial_solve / parallel_solve.max(f64::MIN_POSITIVE)
    );

    if !serial_verified.is_finite() || !parallel_verified.is_finite() {
        return Err("non-finite independently verified residual".into());
    }
    Ok(())
}

examples/fem_structural_pcg_profile.rs (line 557)

fn main() -> Result<(), Box<dyn Error>> {
    let args = Args::parse()?;
    println!(
        "HyBIT {} PCG vector-kernel profile",
        env!("CARGO_PKG_VERSION")
    );
    println!("matrix             : {}", args.matrix.display());
    println!("coordinates        : {}", args.coordinates.display());

    let load_start = Instant::now();
    let (matrix, mm) = read_matrix_market(&args.matrix)?;
    let matrix_load = load_start.elapsed().as_secs_f64();
    let coord_start = Instant::now();
    let coordinates = read_coordinates(&args.coordinates)?;
    let coord_load = coord_start.elapsed().as_secs_f64();
    let matrix_profile = analyze_csr32(&matrix)?;

    println!(
        "Matrix Market      : {:?}, {} input entries -> {} CSR nnz",
        mm.symmetry, mm.input_entries, mm.csr_nnz
    );
    println!(
        "dimensions         : {} x {}",
        matrix_profile.nrows, matrix_profile.ncols
    );
    println!("nnz                : {}", matrix_profile.nnz);
    println!(
        "CSR storage        : {:.3} MiB",
        mib(matrix.storage_bytes())
    );
    println!("matrix load        : {:.3} ms", matrix_load * 1.0e3);
    println!("coordinate nodes   : {}", coordinates.len());
    println!("coordinate load    : {:.3} ms", coord_load * 1.0e3);
    println!("aggregation        : {:?}", args.aggregation);
    println!("target coarse dim  : {}", args.target_coarse_dimension);

    if !matrix_profile.square || !matrix_profile.full_diagonal || !matrix_profile.positive_diagonal
    {
        return Err(
            "structural PCG benchmark requires a square matrix with a complete positive diagonal"
                .into(),
        );
    }
    if matrix.nrows() != coordinates.len() * 3 {
        return Err(format!(
            "matrix/coordinate mismatch: {} matrix rows != {} coordinate nodes * 3",
            matrix.nrows(),
            coordinates.len()
        )
        .into());
    }

    let b = if let Some(path) = args.rhs.as_deref() {
        println!("RHS                : {}", path.display());
        load_rhs(path, matrix.nrows())?
    } else {
        println!("RHS                : generated as b=A*1 (known exact solution)");
        matrix.spmv(&vec![1.0; matrix.ncols()])?
    };

    let aggregate_nodes =
        recommend_rigid_body_aggregate_nodes(coordinates.len(), args.target_coarse_dimension)?;
    let setup_start = Instant::now();
    let preconditioner = match args.aggregation {
        RigidBodyAggregation::Graph => {
            RigidBodyTwoLevelBlockJacobiPreconditioner::from_csr32_graph(
                &matrix,
                &coordinates,
                aggregate_nodes,
            )?
        }
        RigidBodyAggregation::Contiguous => RigidBodyTwoLevelBlockJacobiPreconditioner::from_csr32(
            &matrix,
            &coordinates,
            aggregate_nodes,
        )?,
        RigidBodyAggregation::Auto => unreachable!(),
    };
    let parallel_preconditioner = ParallelRigidBodyTwoLevelPreconditioner::new(&preconditioner)?;
    let parallel_operator = ParallelCsr32Operator::new(&matrix);
    let setup_seconds = setup_start.elapsed().as_secs_f64();

    println!("aggregate target   : {} nodes", aggregate_nodes);
    println!("aggregate count    : {}", preconditioner.aggregate_count());
    println!("coarse dimension   : {}", preconditioner.coarse_dimension());
    println!(
        "prec storage       : {:.3} MiB",
        mib(preconditioner.factor_bytes())
    );
    println!(
        "parallel index     : {:.3} MiB",
        mib(parallel_preconditioner.index_storage_bytes())
    );
    println!("setup              : {:.3} ms", setup_seconds * 1.0e3);
    println!(
        "Rayon threads      : {}",
        parallel_preconditioner.rayon_threads()
    );

    let options = SolverOptions {
        relative_tolerance: args.relative_tolerance,
        absolute_tolerance: 0.0,
        max_iterations: args.max_iterations,
    };

    let mut x_profiled = vec![0.0; matrix.ncols()];
    let (profiled_out, stages, profiled_wall) = profiled_pcg(
        &parallel_operator,
        &parallel_preconditioner,
        &b,
        &mut x_profiled,
        options,
    )?;
    let profiled_verified = verified_relative_residual(&matrix, &b, &x_profiled)?;

    let mut x_control = vec![0.0; matrix.ncols()];
    let mut workspace = PcgWorkspace::new(matrix.nrows());
    let control_start = Instant::now();
    let control_out = pcg_with_workspace(
        &parallel_operator,
        &parallel_preconditioner,
        &b,
        &mut x_control,
        options,
        &mut workspace,
    )?;
    let control_wall = control_start.elapsed().as_secs_f64();
    let control_verified = verified_relative_residual(&matrix, &b, &x_control)?;

    println!();
    println!("Profiled PCG");
    println!("status             : {:?}", profiled_out.status);
    println!("iterations         : {}", profiled_out.iterations);
    println!(
        "reported residual  : {:.6e}",
        profiled_out.final_residual / norm2(&b).max(f64::MIN_POSITIVE)
    );
    println!("verified residual  : {:.6e}", profiled_verified);
    println!("profiled wall      : {:.3} ms", profiled_wall * 1.0e3);
    println!(
        "observed / iter    : {:.3} ms",
        if profiled_out.iterations == 0 {
            0.0
        } else {
            profiled_wall * 1.0e3 / profiled_out.iterations as f64
        }
    );
    println!();
    println!("Actual nested PCG stage timings");
    print_stage(
        "operator A*x",
        stages.operator_seconds,
        stages.operator_calls,
        profiled_wall,
    );
    print_stage(
        "preconditioner",
        stages.preconditioner_seconds,
        stages.preconditioner_calls,
        profiled_wall,
    );
    print_stage(
        "dot reductions",
        stages.dot_seconds,
        stages.dot_calls,
        profiled_wall,
    );
    print_stage(
        "norm reductions",
        stages.norm_seconds,
        stages.norm_calls,
        profiled_wall,
    );
    print_stage(
        "x/r fused update",
        stages.update_x_r_seconds,
        stages.update_x_r_calls,
        profiled_wall,
    );
    print_stage(
        "p update",
        stages.update_p_seconds,
        stages.update_p_calls,
        profiled_wall,
    );
    print_stage(
        "initial residual",
        stages.residual_init_seconds,
        1,
        profiled_wall,
    );
    print_stage("initial p copy", stages.copy_p_seconds, 1, profiled_wall);
    let accounted = stages.accounted_seconds();
    let unaccounted = (profiled_wall - accounted).max(0.0);
    println!(
        "accounted total        : {:9.3} ms  {:5.1}%",
        accounted * 1.0e3,
        100.0 * accounted / profiled_wall.max(f64::MIN_POSITIVE)
    );
    println!(
        "timer/control overhead : {:9.3} ms  {:5.1}%",
        unaccounted * 1.0e3,
        100.0 * unaccounted / profiled_wall.max(f64::MIN_POSITIVE)
    );

    println!();
    println!("Uninstrumented control");
    println!("status             : {:?}", control_out.status);
    println!("iterations         : {}", control_out.iterations);
    println!(
        "reported residual  : {:.6e}",
        control_out.final_residual / norm2(&b).max(f64::MIN_POSITIVE)
    );
    println!("verified residual  : {:.6e}", control_verified);
    println!("solve time         : {:.3} ms", control_wall * 1.0e3);
    println!(
        "profile/control    : {:.3}x",
        profiled_wall / control_wall.max(f64::MIN_POSITIVE)
    );
    println!("note               : profile uses the same serial PCG vector kernels as production; only nested timers are added");

    Ok(())
}

examples/fem_structural_precond_parallel.rs (line 322)

fn main() -> Result<(), Box<dyn Error>> {
    let args = Args::parse()?;
    println!(
        "HyBIT {} serial vs parallel rigid-body preconditioner benchmark",
        env!("CARGO_PKG_VERSION")
    );
    println!("matrix             : {}", args.matrix.display());
    println!("coordinates        : {}", args.coordinates.display());

    let load_start = Instant::now();
    let (matrix, mm) = read_matrix_market(&args.matrix)?;
    let matrix_load = load_start.elapsed().as_secs_f64();
    let coord_start = Instant::now();
    let coordinates = read_coordinates(&args.coordinates)?;
    let coord_load = coord_start.elapsed().as_secs_f64();
    let profile = analyze_csr32(&matrix)?;

    println!(
        "Matrix Market      : {:?}, {} input entries -> {} CSR nnz",
        mm.symmetry, mm.input_entries, mm.csr_nnz
    );
    println!("dimensions         : {} x {}", profile.nrows, profile.ncols);
    println!("nnz                : {}", profile.nnz);
    println!(
        "CSR storage        : {:.3} MiB",
        mib(matrix.storage_bytes())
    );
    println!("matrix load        : {:.3} ms", matrix_load * 1.0e3);
    println!("coordinate nodes   : {}", coordinates.len());
    println!("coordinate load    : {:.3} ms", coord_load * 1.0e3);
    println!("aggregation        : {:?}", args.aggregation);
    println!("target coarse dim  : {}", args.target_coarse_dimension);
    println!("kernel repeats     : {}", args.kernel_repeats);

    if !profile.square || !profile.full_diagonal || !profile.positive_diagonal {
        return Err(
            "structural PCG benchmark requires a square matrix with a complete positive diagonal"
                .into(),
        );
    }
    if matrix.nrows() != coordinates.len() * 3 {
        return Err(format!(
            "matrix/coordinate mismatch: {} matrix rows != {} coordinate nodes * 3",
            matrix.nrows(),
            coordinates.len()
        )
        .into());
    }

    let b = if let Some(path) = args.rhs.as_deref() {
        println!("RHS                : {}", path.display());
        load_rhs(path, matrix.nrows())?
    } else {
        println!("RHS                : generated as b=A*1 (known exact solution)");
        matrix.spmv(&vec![1.0; matrix.ncols()])?
    };

    let aggregate_nodes =
        recommend_rigid_body_aggregate_nodes(coordinates.len(), args.target_coarse_dimension)?;
    let setup_start = Instant::now();
    let preconditioner = match args.aggregation {
        RigidBodyAggregation::Graph => {
            RigidBodyTwoLevelBlockJacobiPreconditioner::from_csr32_graph(
                &matrix,
                &coordinates,
                aggregate_nodes,
            )?
        }
        RigidBodyAggregation::Contiguous => RigidBodyTwoLevelBlockJacobiPreconditioner::from_csr32(
            &matrix,
            &coordinates,
            aggregate_nodes,
        )?,
        RigidBodyAggregation::Auto => unreachable!(),
    };
    let setup_seconds = setup_start.elapsed().as_secs_f64();
    let parallel_preconditioner = ParallelRigidBodyTwoLevelPreconditioner::new(&preconditioner)?;
    let parallel_operator = ParallelCsr32Operator::new(&matrix);

    println!("aggregate target   : {} nodes", aggregate_nodes);
    println!("aggregate count    : {}", preconditioner.aggregate_count());
    println!("coarse dimension   : {}", preconditioner.coarse_dimension());
    println!(
        "prec storage       : {:.3} MiB",
        mib(preconditioner.factor_bytes())
    );
    println!(
        "parallel index     : {:.3} MiB",
        mib(parallel_preconditioner.index_storage_bytes())
    );
    println!("setup              : {:.3} ms", setup_seconds * 1.0e3);
    println!(
        "Rayon threads      : {}",
        parallel_preconditioner.rayon_threads()
    );

    let mut zs = vec![0.0; matrix.nrows()];
    let mut zp = vec![0.0; matrix.nrows()];
    preconditioner.apply(&b, &mut zs)?;
    parallel_preconditioner.apply(&b, &mut zp)?;
    let scale = zs.iter().fold(1.0f64, |m, &v| m.max(v.abs()));
    let max_diff = zs
        .iter()
        .zip(&zp)
        .fold(0.0f64, |m, (&a, &b)| m.max((a - b).abs()));
    if max_diff > 1.0e-11 * scale {
        return Err(format!(
            "serial/parallel preconditioner mismatch: max diff={max_diff:e}, scale={scale:e}"
        )
        .into());
    }

    let serial_prec_start = Instant::now();
    for _ in 0..args.kernel_repeats {
        preconditioner.apply(&b, &mut zs)?;
    }
    let serial_prec = serial_prec_start.elapsed().as_secs_f64();

    let parallel_prec_start = Instant::now();
    for _ in 0..args.kernel_repeats {
        parallel_preconditioner.apply(&b, &mut zp)?;
    }
    let parallel_prec = parallel_prec_start.elapsed().as_secs_f64();

    println!();
    println!("Preconditioner microbenchmark");
    println!(
        "serial / apply     : {:.3} ms",
        ms_per(serial_prec, args.kernel_repeats)
    );
    println!(
        "parallel / apply   : {:.3} ms",
        ms_per(parallel_prec, args.kernel_repeats)
    );
    println!(
        "precond speedup    : {:.3}x",
        serial_prec / parallel_prec.max(f64::MIN_POSITIVE)
    );

    let options = SolverOptions {
        relative_tolerance: args.relative_tolerance,
        absolute_tolerance: 0.0,
        max_iterations: args.max_iterations,
    };

    let mut xs = vec![0.0; matrix.ncols()];
    let mut ws = PcgWorkspace::new(matrix.nrows());
    let serial_solve_start = Instant::now();
    let serial_out = pcg_with_workspace(
        &parallel_operator,
        &preconditioner,
        &b,
        &mut xs,
        options,
        &mut ws,
    )?;
    let serial_solve = serial_solve_start.elapsed().as_secs_f64();
    let serial_verified = verified_relative_residual(&matrix, &b, &xs)?;

    let mut xp = vec![0.0; matrix.ncols()];
    let mut wp = PcgWorkspace::new(matrix.nrows());
    let parallel_solve_start = Instant::now();
    let parallel_out = pcg_with_workspace(
        &parallel_operator,
        &parallel_preconditioner,
        &b,
        &mut xp,
        options,
        &mut wp,
    )?;
    let parallel_solve = parallel_solve_start.elapsed().as_secs_f64();
    let parallel_verified = verified_relative_residual(&matrix, &b, &xp)?;

    println!();
    println!("[1/2] Parallel CSR + serial preconditioner");
    println!("status             : {:?}", serial_out.status);
    println!("iterations         : {}", serial_out.iterations);
    println!(
        "reported residual  : {:.6e}",
        serial_out.final_residual / norm2(&b).max(f64::MIN_POSITIVE)
    );
    println!("verified residual  : {:.6e}", serial_verified);
    println!("solve time         : {:.3} ms", serial_solve * 1.0e3);

    println!();
    println!("[2/2] Parallel CSR + parallel preconditioner");
    println!("status             : {:?}", parallel_out.status);
    println!("iterations         : {}", parallel_out.iterations);
    println!(
        "reported residual  : {:.6e}",
        parallel_out.final_residual / norm2(&b).max(f64::MIN_POSITIVE)
    );
    println!("verified residual  : {:.6e}", parallel_verified);
    println!("solve time         : {:.3} ms", parallel_solve * 1.0e3);

    println!();
    println!("Comparison");
    println!(
        "iteration ratio    : {:.3} (parallel prec / serial prec)",
        parallel_out.iterations as f64 / serial_out.iterations.max(1) as f64
    );
    println!(
        "solve-time ratio   : {:.3} (parallel prec / serial prec)",
        parallel_solve / serial_solve.max(f64::MIN_POSITIVE)
    );
    println!(
        "PCG speedup        : {:.3}x",
        serial_solve / parallel_solve.max(f64::MIN_POSITIVE)
    );

    if !serial_verified.is_finite() || !parallel_verified.is_finite() {
        return Err("non-finite independently verified residual".into());
    }
    Ok(())
}

examples/fem_structural_prepared.rs (line 402)

fn main() -> Result<(), Box<dyn Error>> {
    let args = Args::parse()?;
    println!(
        "HyBIT {} prepared structural solve-many FEM benchmark",
        env!("CARGO_PKG_VERSION")
    );
    println!("matrix             : {}", args.matrix.display());
    println!("coordinates        : {}", args.coordinates.display());

    let load_start = Instant::now();
    let (matrix, mm) = read_matrix_market(&args.matrix)?;
    let matrix_load_seconds = load_start.elapsed().as_secs_f64();
    let coord_start = Instant::now();
    let coordinates = read_coordinates(&args.coordinates)?;
    let coord_load_seconds = coord_start.elapsed().as_secs_f64();
    let profile = analyze_csr32(&matrix)?;

    println!(
        "Matrix Market      : {:?}, {} input entries -> {} CSR nnz",
        mm.symmetry, mm.input_entries, mm.csr_nnz
    );
    println!("dimensions         : {} x {}", profile.nrows, profile.ncols);
    println!("nnz                : {}", profile.nnz);
    println!(
        "CSR storage        : {:.3} MiB",
        mib(matrix.storage_bytes())
    );
    println!("matrix load        : {:.3} ms", matrix_load_seconds * 1.0e3);
    println!("coordinate nodes   : {}", coordinates.len());
    println!("coordinate load    : {:.3} ms", coord_load_seconds * 1.0e3);

    if !profile.square || !profile.full_diagonal || !profile.positive_diagonal {
        return Err(
            "prepared structural PCG requires a square matrix with a complete positive diagonal"
                .into(),
        );
    }
    if matrix.nrows() != coordinates.len() * 3 {
        return Err(format!(
            "matrix/coordinate mismatch: {} matrix rows != {} coordinate nodes * 3",
            matrix.nrows(),
            coordinates.len()
        )
        .into());
    }

    let generated_rhs = args.rhs.is_none();
    let b = if let Some(path) = args.rhs.as_deref() {
        println!("RHS                : {}", path.display());
        load_rhs(path, matrix.nrows())?
    } else {
        println!("RHS                : generated as b=A*1 (known exact solution)");
        let ones = vec![1.0; matrix.ncols()];
        matrix.spmv(&ones)?
    };
    println!(
        "repeats            : {} (fresh zero initial guess each solve)",
        args.repeats
    );

    let mut solver = HybitSolver::new();
    solver.set_options(SolverOptions {
        relative_tolerance: args.relative_tolerance,
        absolute_tolerance: 0.0,
        max_iterations: args.max_iterations,
    })?;
    solver.set_structural_options(StructuralOptions {
        target_coarse_dimension: args.target_coarse_dimension,
        aggregation: args.aggregation,
        spmv_policy: args.spmv_policy,
        preconditioner_policy: args.preconditioner_policy,
        pcg_vector_policy: args.pcg_vector_policy,
    })?;

    let analysis = solver.analyze_csr32(&matrix)?;
    let mut prepared = solver.prepare_structural_csr32(&matrix, &analysis, &coordinates)?;

    println!("policy             : StructuralAuto/RigidBodyTwoLevel/Prepared");
    println!("aggregation        : {:?}", prepared.aggregation());
    println!("SpMV policy        : {:?}", prepared.spmv_policy());
    println!(
        "precond policy     : {:?}",
        prepared.structural_preconditioner_policy()
    );
    println!("PCG vector policy  : {:?}", prepared.pcg_vector_policy());
    if prepared.parallel_spmv_enabled()
        || prepared.parallel_preconditioner_enabled()
        || prepared.parallel_pcg_vectors_enabled()
    {
        println!(
            "Rayon threads      : {}",
            ParallelCsr32Operator::new(&matrix).rayon_threads()
        );
    }
    if prepared.parallel_preconditioner_enabled() {
        println!(
            "parallel index     : {:.3} MiB",
            mib(prepared.parallel_preconditioner_index_bytes())
        );
    }
    println!("target coarse dim  : {}", args.target_coarse_dimension);
    println!(
        "aggregate nodes    : {} (auto-selected)",
        prepared.aggregate_nodes()
    );
    println!("aggregate count    : {}", prepared.aggregate_count());
    println!(
        "aggregate min/max  : {} / {} nodes",
        prepared.min_aggregate_nodes(),
        prepared.max_aggregate_nodes()
    );
    println!("modes/aggregate    : 6");
    println!("coarse dimension   : {}", prepared.coarse_dimension());
    println!(
        "base factor        : {:.3} MiB",
        mib(prepared.base_factor_bytes())
    );
    println!(
        "coarse factor      : {:.3} MiB",
        mib(prepared.coarse_factor_bytes())
    );
    println!(
        "geometry storage   : {:.3} MiB",
        mib(prepared.geometry_bytes())
    );
    println!(
        "total prec storage : {:.3} MiB",
        mib(prepared.preconditioner_bytes())
    );
    println!(
        "Krylov workspace   : {:.3} MiB",
        mib(prepared.krylov_workspace_bytes())
    );
    println!(
        "analysis once      : {:.3} ms",
        prepared.analysis_seconds() * 1.0e3
    );
    println!(
        "prepare once       : {:.3} ms",
        prepared.prepare_seconds() * 1.0e3
    );

    let mut total_solve_seconds = 0.0f64;
    for repetition in 1..=args.repeats {
        let mut x = vec![0.0; matrix.ncols()];
        let wall_start = Instant::now();
        let report = prepared.solve(&matrix, &b, &mut x)?;
        let wall_seconds = wall_start.elapsed().as_secs_f64();
        let verified = verified_relative_residual(&matrix, &b, &x)?;
        total_solve_seconds += report.solve_seconds;

        println!();
        println!("Solve #{repetition}");
        println!("status             : {:?}", report.status);
        println!("preconditioner     : {:?}", report.preconditioner);
        println!("reused             : {}", report.preconditioner_reused);
        println!("sequence           : {}", report.solve_sequence);
        println!("iterations         : {}", report.iterations);
        println!("reported residual  : {:.6e}", report.relative_residual);
        println!("verified residual  : {:.6e}", verified);
        if generated_rhs {
            println!("relative x error   : {:.6e}", relative_error_to_ones(&x));
        }
        println!(
            "analysis charged   : {:.3} ms",
            report.analysis_seconds * 1.0e3
        );
        println!(
            "prepare charged    : {:.3} ms",
            report.prepare_seconds * 1.0e3
        );
        println!(
            "solve time         : {:.3} ms",
            report.solve_seconds * 1.0e3
        );
        println!("call wall          : {:.3} ms", wall_seconds * 1.0e3);

        if repetition == 1 && report.preconditioner_reused {
            return Err("first prepared structural solve unexpectedly reported reuse".into());
        }
        if repetition > 1 && !report.preconditioner_reused {
            return Err(
                "subsequent prepared structural solve did not report preconditioner reuse".into(),
            );
        }
        if repetition > 1 && (report.analysis_seconds != 0.0 || report.prepare_seconds != 0.0) {
            return Err(
                "subsequent prepared structural solve was charged reusable setup cost".into(),
            );
        }
        if !verified.is_finite() {
            return Err("non-finite independently verified residual".into());
        }
    }

    println!();
    println!("Prepared reuse summary");
    println!("solve count        : {}", prepared.solve_count());
    println!(
        "one-time setup     : {:.3} ms",
        (prepared.analysis_seconds() + prepared.prepare_seconds()) * 1.0e3
    );
    println!("sum solve time     : {:.3} ms", total_solve_seconds * 1.0e3);
    println!(
        "amortized total   : {:.3} ms/solve",
        ((prepared.analysis_seconds() + prepared.prepare_seconds() + total_solve_seconds)
            / args.repeats as f64)
            * 1.0e3
    );

    Ok(())
}

pub fn matrix(&self) -> &'a Csr32Matrix

pub fn rayon_threads(&self) -> usize

Examples found in repository

examples/fem_structural_auto.rs (line 389)

fn main() -> Result<(), Box<dyn Error>> {
    let args = Args::parse()?;
    println!(
        "HyBIT {} structural-auto FEM benchmark",
        env!("CARGO_PKG_VERSION")
    );
    println!("matrix             : {}", args.matrix.display());
    println!("coordinates        : {}", args.coordinates.display());

    let load_start = Instant::now();
    let (matrix, mm) = read_matrix_market(&args.matrix)?;
    let matrix_load_seconds = load_start.elapsed().as_secs_f64();
    let coord_start = Instant::now();
    let coordinates = read_coordinates(&args.coordinates)?;
    let coord_load_seconds = coord_start.elapsed().as_secs_f64();
    let profile = analyze_csr32(&matrix)?;

    println!(
        "Matrix Market      : {:?}, {} input entries -> {} CSR nnz",
        mm.symmetry, mm.input_entries, mm.csr_nnz
    );
    println!("dimensions         : {} x {}", profile.nrows, profile.ncols);
    println!("nnz                : {}", profile.nnz);
    println!(
        "CSR storage        : {:.3} MiB",
        mib(matrix.storage_bytes())
    );
    println!("matrix load        : {:.3} ms", matrix_load_seconds * 1.0e3);
    println!("coordinate nodes   : {}", coordinates.len());
    println!("coordinate load    : {:.3} ms", coord_load_seconds * 1.0e3);

    if !profile.square || !profile.full_diagonal || !profile.positive_diagonal {
        return Err(
            "rigid-body two-level PCG requires a square matrix with a complete positive diagonal"
                .into(),
        );
    }
    if matrix.nrows() != coordinates.len() * 3 {
        return Err(format!(
            "matrix/coordinate mismatch: {} matrix rows != {} coordinate nodes * 3",
            matrix.nrows(),
            coordinates.len()
        )
        .into());
    }

    let generated_rhs = args.rhs.is_none();
    let b = if let Some(path) = args.rhs.as_deref() {
        println!("RHS                : {}", path.display());
        load_rhs(path, matrix.nrows())?
    } else {
        println!("RHS                : generated as b=A*1 (known exact solution)");
        let ones = vec![1.0; matrix.ncols()];
        matrix.spmv(&ones)?
    };
    let mut solver = HybitSolver::new();
    solver.set_options(SolverOptions {
        relative_tolerance: args.relative_tolerance,
        absolute_tolerance: 0.0,
        max_iterations: args.max_iterations,
    })?;
    solver.set_structural_options(StructuralOptions {
        target_coarse_dimension: args.target_coarse_dimension,
        aggregation: args.aggregation,
        spmv_policy: args.spmv_policy,
        preconditioner_policy: args.preconditioner_policy,
        pcg_vector_policy: args.pcg_vector_policy,
    })?;

    let analysis = solver.analyze_csr32(&matrix)?;
    let mut prepared = solver.prepare_structural_csr32(&matrix, &analysis, &coordinates)?;
    println!("policy             : StructuralAuto/RigidBodyTwoLevel");
    println!("aggregation        : {:?}", prepared.aggregation());
    println!("SpMV policy        : {:?}", prepared.spmv_policy());
    println!(
        "precond policy     : {:?}",
        prepared.structural_preconditioner_policy()
    );
    println!("PCG vector policy  : {:?}", prepared.pcg_vector_policy());
    if prepared.parallel_spmv_enabled()
        || prepared.parallel_preconditioner_enabled()
        || prepared.parallel_pcg_vectors_enabled()
    {
        println!(
            "Rayon threads      : {}",
            ParallelCsr32Operator::new(&matrix).rayon_threads()
        );
    }
    if prepared.parallel_preconditioner_enabled() {
        println!(
            "parallel index     : {:.3} MiB",
            mib(prepared.parallel_preconditioner_index_bytes())
        );
    }
    println!("target coarse dim  : {}", args.target_coarse_dimension);
    println!(
        "aggregate nodes    : {} (auto-selected)",
        prepared.aggregate_nodes()
    );
    println!("fine block size    : 3");
    println!("aggregate count    : {}", prepared.aggregate_count());
    println!(
        "aggregate min/max  : {} / {} nodes",
        prepared.min_aggregate_nodes(),
        prepared.max_aggregate_nodes()
    );
    println!("modes/aggregate    : 6");
    println!("coarse dimension   : {}", prepared.coarse_dimension());
    println!(
        "base factor        : {:.3} MiB",
        mib(prepared.base_factor_bytes())
    );
    println!(
        "coarse factor      : {:.3} MiB",
        mib(prepared.coarse_factor_bytes())
    );
    println!(
        "geometry storage   : {:.3} MiB",
        mib(prepared.geometry_bytes())
    );
    println!(
        "total prec storage : {:.3} MiB",
        mib(prepared.preconditioner_bytes())
    );
    println!(
        "analysis           : {:.3} ms",
        prepared.analysis_seconds() * 1.0e3
    );
    println!(
        "prepare            : {:.3} ms",
        prepared.prepare_seconds() * 1.0e3
    );

    let mut x = vec![0.0; matrix.ncols()];
    let report = prepared.solve(&matrix, &b, &mut x)?;
    let verified = verified_relative_residual(&matrix, &b, &x)?;

    println!();
    println!("Structural Auto API: 3x3 Block-Jacobi + six rigid-body coarse modes");
    println!("status             : {:?}", report.status);
    println!("preconditioner     : {:?}", report.preconditioner);
    println!("iterations         : {}", report.iterations);
    println!("reported residual  : {:.6e}", report.relative_residual);
    println!("verified residual  : {:.6e}", verified);
    if generated_rhs {
        println!("relative x error   : {:.6e}", relative_error_to_ones(&x));
    }
    println!(
        "solve time         : {:.3} ms",
        report.solve_seconds * 1.0e3
    );
    println!(
        "total setup+solve  : {:.3} ms",
        (report.setup_seconds + report.solve_seconds) * 1.0e3
    );

    if !verified.is_finite() {
        return Err("non-finite independently verified residual".into());
    }
    Ok(())
}

examples/fem_structural_spmv.rs (line 329)

fn main() -> Result<(), Box<dyn Error>> {
    let args = Args::parse()?;
    println!(
        "HyBIT {} serial vs parallel CSR structural benchmark",
        env!("CARGO_PKG_VERSION")
    );
    println!("matrix             : {}", args.matrix.display());
    println!("coordinates        : {}", args.coordinates.display());

    let load_start = Instant::now();
    let (matrix, mm) = read_matrix_market(&args.matrix)?;
    let matrix_load = load_start.elapsed().as_secs_f64();
    let coord_start = Instant::now();
    let coordinates = read_coordinates(&args.coordinates)?;
    let coord_load = coord_start.elapsed().as_secs_f64();
    let profile = analyze_csr32(&matrix)?;

    println!(
        "Matrix Market      : {:?}, {} input entries -> {} CSR nnz",
        mm.symmetry, mm.input_entries, mm.csr_nnz
    );
    println!("dimensions         : {} x {}", profile.nrows, profile.ncols);
    println!("nnz                : {}", profile.nnz);
    println!(
        "CSR storage        : {:.3} MiB",
        mib(matrix.storage_bytes())
    );
    println!("matrix load        : {:.3} ms", matrix_load * 1.0e3);
    println!("coordinate nodes   : {}", coordinates.len());
    println!("coordinate load    : {:.3} ms", coord_load * 1.0e3);
    println!("aggregation        : {:?}", args.aggregation);
    println!("target coarse dim  : {}", args.target_coarse_dimension);
    println!("kernel repeats     : {}", args.kernel_repeats);

    if !profile.square || !profile.full_diagonal || !profile.positive_diagonal {
        return Err(
            "structural PCG benchmark requires a square matrix with a complete positive diagonal"
                .into(),
        );
    }
    if matrix.nrows() != coordinates.len() * 3 {
        return Err(format!(
            "matrix/coordinate mismatch: {} matrix rows != {} coordinate nodes * 3",
            matrix.nrows(),
            coordinates.len()
        )
        .into());
    }

    let b = if let Some(path) = args.rhs.as_deref() {
        println!("RHS                : {}", path.display());
        load_rhs(path, matrix.nrows())?
    } else {
        println!("RHS                : generated as b=A*1 (known exact solution)");
        matrix.spmv(&vec![1.0; matrix.ncols()])?
    };

    let aggregate_nodes =
        recommend_rigid_body_aggregate_nodes(coordinates.len(), args.target_coarse_dimension)?;
    let setup_start = Instant::now();
    let preconditioner = match args.aggregation {
        RigidBodyAggregation::Graph => {
            RigidBodyTwoLevelBlockJacobiPreconditioner::from_csr32_graph(
                &matrix,
                &coordinates,
                aggregate_nodes,
            )?
        }
        RigidBodyAggregation::Contiguous => RigidBodyTwoLevelBlockJacobiPreconditioner::from_csr32(
            &matrix,
            &coordinates,
            aggregate_nodes,
        )?,
        RigidBodyAggregation::Auto => unreachable!(),
    };
    let setup_seconds = setup_start.elapsed().as_secs_f64();
    let parallel = ParallelCsr32Operator::new(&matrix);

    println!("aggregate target   : {} nodes", aggregate_nodes);
    println!("aggregate count    : {}", preconditioner.aggregate_count());
    println!("coarse dimension   : {}", preconditioner.coarse_dimension());
    println!(
        "prec storage       : {:.3} MiB",
        mib(preconditioner.factor_bytes())
    );
    println!("setup              : {:.3} ms", setup_seconds * 1.0e3);
    println!("Rayon threads      : {}", parallel.rayon_threads());

    let mut ys = vec![0.0; matrix.nrows()];
    let mut yp = vec![0.0; matrix.nrows()];
    matrix.apply(&b, &mut ys)?;
    parallel.apply(&b, &mut yp)?;
    if ys != yp {
        return Err("serial and parallel CSR kernels produced different results".into());
    }

    let serial_start = Instant::now();
    for _ in 0..args.kernel_repeats {
        matrix.apply(&b, &mut ys)?;
    }
    let serial_spmv = serial_start.elapsed().as_secs_f64();

    let parallel_start = Instant::now();
    for _ in 0..args.kernel_repeats {
        parallel.apply(&b, &mut yp)?;
    }
    let parallel_spmv = parallel_start.elapsed().as_secs_f64();

    println!();
    println!("CSR SpMV microbenchmark");
    println!(
        "serial / call      : {:.3} ms",
        ms_per(serial_spmv, args.kernel_repeats)
    );
    println!(
        "parallel / call    : {:.3} ms",
        ms_per(parallel_spmv, args.kernel_repeats)
    );
    println!(
        "SpMV speedup       : {:.3}x",
        serial_spmv / parallel_spmv.max(f64::MIN_POSITIVE)
    );

    let options = SolverOptions {
        relative_tolerance: args.relative_tolerance,
        absolute_tolerance: 0.0,
        max_iterations: args.max_iterations,
    };

    let mut xs = vec![0.0; matrix.ncols()];
    let mut ws = PcgWorkspace::new(matrix.nrows());
    let serial_solve_start = Instant::now();
    let serial_out = pcg_with_workspace(&matrix, &preconditioner, &b, &mut xs, options, &mut ws)?;
    let serial_solve = serial_solve_start.elapsed().as_secs_f64();
    let serial_verified = verified_relative_residual(&matrix, &b, &xs)?;

    let mut xp = vec![0.0; matrix.ncols()];
    let mut wp = PcgWorkspace::new(matrix.nrows());
    let parallel_solve_start = Instant::now();
    let parallel_out =
        pcg_with_workspace(&parallel, &preconditioner, &b, &mut xp, options, &mut wp)?;
    let parallel_solve = parallel_solve_start.elapsed().as_secs_f64();
    let parallel_verified = verified_relative_residual(&matrix, &b, &xp)?;

    println!();
    println!("[1/2] Serial CSR PCG");
    println!("status             : {:?}", serial_out.status);
    println!("iterations         : {}", serial_out.iterations);
    println!(
        "reported residual  : {:.6e}",
        serial_out.final_residual / norm2(&b).max(f64::MIN_POSITIVE)
    );
    println!("verified residual  : {:.6e}", serial_verified);
    println!("solve time         : {:.3} ms", serial_solve * 1.0e3);

    println!();
    println!("[2/2] Parallel CSR PCG");
    println!("status             : {:?}", parallel_out.status);
    println!("iterations         : {}", parallel_out.iterations);
    println!(
        "reported residual  : {:.6e}",
        parallel_out.final_residual / norm2(&b).max(f64::MIN_POSITIVE)
    );
    println!("verified residual  : {:.6e}", parallel_verified);
    println!("solve time         : {:.3} ms", parallel_solve * 1.0e3);

    println!();
    println!("Comparison");
    println!(
        "iteration ratio    : {:.3} (parallel / serial)",
        parallel_out.iterations as f64 / serial_out.iterations.max(1) as f64
    );
    println!(
        "solve-time ratio   : {:.3} (parallel / serial)",
        parallel_solve / serial_solve.max(f64::MIN_POSITIVE)
    );
    println!(
        "PCG speedup        : {:.3}x",
        serial_solve / parallel_solve.max(f64::MIN_POSITIVE)
    );

    if !serial_verified.is_finite() || !parallel_verified.is_finite() {
        return Err("non-finite independently verified residual".into());
    }
    Ok(())
}

examples/fem_structural_prepared.rs (line 402)

fn main() -> Result<(), Box<dyn Error>> {
    let args = Args::parse()?;
    println!(
        "HyBIT {} prepared structural solve-many FEM benchmark",
        env!("CARGO_PKG_VERSION")
    );
    println!("matrix             : {}", args.matrix.display());
    println!("coordinates        : {}", args.coordinates.display());

    let load_start = Instant::now();
    let (matrix, mm) = read_matrix_market(&args.matrix)?;
    let matrix_load_seconds = load_start.elapsed().as_secs_f64();
    let coord_start = Instant::now();
    let coordinates = read_coordinates(&args.coordinates)?;
    let coord_load_seconds = coord_start.elapsed().as_secs_f64();
    let profile = analyze_csr32(&matrix)?;

    println!(
        "Matrix Market      : {:?}, {} input entries -> {} CSR nnz",
        mm.symmetry, mm.input_entries, mm.csr_nnz
    );
    println!("dimensions         : {} x {}", profile.nrows, profile.ncols);
    println!("nnz                : {}", profile.nnz);
    println!(
        "CSR storage        : {:.3} MiB",
        mib(matrix.storage_bytes())
    );
    println!("matrix load        : {:.3} ms", matrix_load_seconds * 1.0e3);
    println!("coordinate nodes   : {}", coordinates.len());
    println!("coordinate load    : {:.3} ms", coord_load_seconds * 1.0e3);

    if !profile.square || !profile.full_diagonal || !profile.positive_diagonal {
        return Err(
            "prepared structural PCG requires a square matrix with a complete positive diagonal"
                .into(),
        );
    }
    if matrix.nrows() != coordinates.len() * 3 {
        return Err(format!(
            "matrix/coordinate mismatch: {} matrix rows != {} coordinate nodes * 3",
            matrix.nrows(),
            coordinates.len()
        )
        .into());
    }

    let generated_rhs = args.rhs.is_none();
    let b = if let Some(path) = args.rhs.as_deref() {
        println!("RHS                : {}", path.display());
        load_rhs(path, matrix.nrows())?
    } else {
        println!("RHS                : generated as b=A*1 (known exact solution)");
        let ones = vec![1.0; matrix.ncols()];
        matrix.spmv(&ones)?
    };
    println!(
        "repeats            : {} (fresh zero initial guess each solve)",
        args.repeats
    );

    let mut solver = HybitSolver::new();
    solver.set_options(SolverOptions {
        relative_tolerance: args.relative_tolerance,
        absolute_tolerance: 0.0,
        max_iterations: args.max_iterations,
    })?;
    solver.set_structural_options(StructuralOptions {
        target_coarse_dimension: args.target_coarse_dimension,
        aggregation: args.aggregation,
        spmv_policy: args.spmv_policy,
        preconditioner_policy: args.preconditioner_policy,
        pcg_vector_policy: args.pcg_vector_policy,
    })?;

    let analysis = solver.analyze_csr32(&matrix)?;
    let mut prepared = solver.prepare_structural_csr32(&matrix, &analysis, &coordinates)?;

    println!("policy             : StructuralAuto/RigidBodyTwoLevel/Prepared");
    println!("aggregation        : {:?}", prepared.aggregation());
    println!("SpMV policy        : {:?}", prepared.spmv_policy());
    println!(
        "precond policy     : {:?}",
        prepared.structural_preconditioner_policy()
    );
    println!("PCG vector policy  : {:?}", prepared.pcg_vector_policy());
    if prepared.parallel_spmv_enabled()
        || prepared.parallel_preconditioner_enabled()
        || prepared.parallel_pcg_vectors_enabled()
    {
        println!(
            "Rayon threads      : {}",
            ParallelCsr32Operator::new(&matrix).rayon_threads()
        );
    }
    if prepared.parallel_preconditioner_enabled() {
        println!(
            "parallel index     : {:.3} MiB",
            mib(prepared.parallel_preconditioner_index_bytes())
        );
    }
    println!("target coarse dim  : {}", args.target_coarse_dimension);
    println!(
        "aggregate nodes    : {} (auto-selected)",
        prepared.aggregate_nodes()
    );
    println!("aggregate count    : {}", prepared.aggregate_count());
    println!(
        "aggregate min/max  : {} / {} nodes",
        prepared.min_aggregate_nodes(),
        prepared.max_aggregate_nodes()
    );
    println!("modes/aggregate    : 6");
    println!("coarse dimension   : {}", prepared.coarse_dimension());
    println!(
        "base factor        : {:.3} MiB",
        mib(prepared.base_factor_bytes())
    );
    println!(
        "coarse factor      : {:.3} MiB",
        mib(prepared.coarse_factor_bytes())
    );
    println!(
        "geometry storage   : {:.3} MiB",
        mib(prepared.geometry_bytes())
    );
    println!(
        "total prec storage : {:.3} MiB",
        mib(prepared.preconditioner_bytes())
    );
    println!(
        "Krylov workspace   : {:.3} MiB",
        mib(prepared.krylov_workspace_bytes())
    );
    println!(
        "analysis once      : {:.3} ms",
        prepared.analysis_seconds() * 1.0e3
    );
    println!(
        "prepare once       : {:.3} ms",
        prepared.prepare_seconds() * 1.0e3
    );

    let mut total_solve_seconds = 0.0f64;
    for repetition in 1..=args.repeats {
        let mut x = vec![0.0; matrix.ncols()];
        let wall_start = Instant::now();
        let report = prepared.solve(&matrix, &b, &mut x)?;
        let wall_seconds = wall_start.elapsed().as_secs_f64();
        let verified = verified_relative_residual(&matrix, &b, &x)?;
        total_solve_seconds += report.solve_seconds;

        println!();
        println!("Solve #{repetition}");
        println!("status             : {:?}", report.status);
        println!("preconditioner     : {:?}", report.preconditioner);
        println!("reused             : {}", report.preconditioner_reused);
        println!("sequence           : {}", report.solve_sequence);
        println!("iterations         : {}", report.iterations);
        println!("reported residual  : {:.6e}", report.relative_residual);
        println!("verified residual  : {:.6e}", verified);
        if generated_rhs {
            println!("relative x error   : {:.6e}", relative_error_to_ones(&x));
        }
        println!(
            "analysis charged   : {:.3} ms",
            report.analysis_seconds * 1.0e3
        );
        println!(
            "prepare charged    : {:.3} ms",
            report.prepare_seconds * 1.0e3
        );
        println!(
            "solve time         : {:.3} ms",
            report.solve_seconds * 1.0e3
        );
        println!("call wall          : {:.3} ms", wall_seconds * 1.0e3);

        if repetition == 1 && report.preconditioner_reused {
            return Err("first prepared structural solve unexpectedly reported reuse".into());
        }
        if repetition > 1 && !report.preconditioner_reused {
            return Err(
                "subsequent prepared structural solve did not report preconditioner reuse".into(),
            );
        }
        if repetition > 1 && (report.analysis_seconds != 0.0 || report.prepare_seconds != 0.0) {
            return Err(
                "subsequent prepared structural solve was charged reusable setup cost".into(),
            );
        }
        if !verified.is_finite() {
            return Err("non-finite independently verified residual".into());
        }
    }

    println!();
    println!("Prepared reuse summary");
    println!("solve count        : {}", prepared.solve_count());
    println!(
        "one-time setup     : {:.3} ms",
        (prepared.analysis_seconds() + prepared.prepare_seconds()) * 1.0e3
    );
    println!("sum solve time     : {:.3} ms", total_solve_seconds * 1.0e3);
    println!(
        "amortized total   : {:.3} ms/solve",
        ((prepared.analysis_seconds() + prepared.prepare_seconds() + total_solve_seconds)
            / args.repeats as f64)
            * 1.0e3
    );

    Ok(())
}

Trait Implementations

impl<'a> Clone for ParallelCsr32Operator<'a>

fn clone(&self) -> ParallelCsr32Operator<'a>

Returns a duplicate of the value.

fn clone_from(&mut self, source: &Self)

Performs copy-assignment from source.

impl<'a> Copy for ParallelCsr32Operator<'a>

impl<'a> Debug for ParallelCsr32Operator<'a>

fn fmt(&self, f: &mut Formatter<'_>) -> Result<(), Error>

Formats the value using the given formatter.

impl LinearOperator for ParallelCsr32Operator<'_>

fn rows(&self) -> usize

fn cols(&self) -> usize

fn apply(&self, x: &[f64], y: &mut [f64]) -> Result<(), HybitError>