PETSc - MatMultScale？矩阵 X 向量 X 标量

Question

我正在使用PETSc，我想做一些类似的事情，

我知道我能做到:

Mat A
Vec x,y

MatMult(A,x,y)
VecScale(y,0.5)

我只是好奇是否有一个函数可以一次完成所有这些操作。看起来它会节省一个循环。

MatMultScale(A,x,0.5,y)

有这样的功能吗？

Answer 1

此函数(或任何接近的函数)似乎不在 functions operating on Mat 的列表中。所以对你的问题的简短回答是……不。

如果您经常使用 $y=\frac12 Ax$，解决方案是使用 MatScale(A,0.5); 一次性缩放矩阵。

这样的功能有用吗？检查这一点的一种方法是使用 petsc 的 -log_summary 选项来获取一些分析信息。如果您的矩阵很密集，您会发现 MatMult() 花费的时间比 VecScale() 花费的时间长得多。仅当处理稀疏矩阵且每行有一些非空项时，这个问题才有意义。

这是一个测试它的代码，使用 2xIdentity 作为矩阵:

static char help[] = "Tests solving linear system on 0 by 0 matrix.\n\n";

#include <petscksp.h>

#undef __FUNCT__
#define __FUNCT__ "main"
int main(int argc,char **args)
{
    Vec            x, y;      
    Mat            A;           
    PetscReal      alpha=0.5;  
    PetscErrorCode ierr;
    PetscInt n=42;

    PetscInitialize(&argc,&args,(char*)0,help);
    ierr = PetscOptionsGetInt(NULL,"-n",&n,NULL);CHKERRQ(ierr);

    /* Create the vector*/
    ierr = VecCreate(PETSC_COMM_WORLD,&x);CHKERRQ(ierr);
    ierr = VecSetSizes(x,PETSC_DECIDE,n);CHKERRQ(ierr);
    ierr = VecSetFromOptions(x);CHKERRQ(ierr);
    ierr = VecDuplicate(x,&y);CHKERRQ(ierr);

    /*
     Create matrix.  When using MatCreate(), the matrix format can
     be specified at runtime.

     Performance tuning note:  For problems of substantial size,
     preallocation of matrix memory is crucial for attaining good
     performance. See the matrix chapter of the users manual for details.
     */
    ierr = MatCreate(PETSC_COMM_WORLD,&A);CHKERRQ(ierr);
    ierr = MatSetSizes(A,PETSC_DECIDE,PETSC_DECIDE,n,n);CHKERRQ(ierr);
    ierr = MatSetFromOptions(A);CHKERRQ(ierr);
    ierr = MatSetUp(A);CHKERRQ(ierr);

    /*
This matrix is diagonal, two times identity
should have preallocated, shame
     */
    PetscInt i,col;
    PetscScalar value=2.0;
    for (i=0; i<n; i++) {
        col=i;
        ierr   = MatSetValues(A,1,&i,1,&col,&value,INSERT_VALUES);CHKERRQ(ierr);
    }

    ierr = MatAssemblyBegin(A,MAT_FINAL_ASSEMBLY);CHKERRQ(ierr);
    ierr = MatAssemblyEnd(A,MAT_FINAL_ASSEMBLY);CHKERRQ(ierr);

    /*
     let's do this 42 times for nothing :
     */
    for(i=0;i<42;i++){
        ierr = MatMult(A,x,y);CHKERRQ(ierr);
        ierr = VecScale(y,alpha);CHKERRQ(ierr);
    }

    ierr = VecDestroy(&x);CHKERRQ(ierr);
    ierr = VecDestroy(&y);CHKERRQ(ierr); 
    ierr = MatDestroy(&A);CHKERRQ(ierr);

    ierr = PetscFinalize();
    return 0;
}

生成文件:

include ${PETSC_DIR}/conf/variables
include ${PETSC_DIR}/conf/rules
include ${PETSC_DIR}/conf/test

CLINKER=g++

all : ex1

ex1 : main.o chkopts
    ${CLINKER} -w -o main main.o ${PETSC_LIB}
    ${RM} main.o

run :
    mpirun -np 2 main -n 10000000 -log_summary -help -mat_type mpiaij

这里生成的两行-log_summary可以回答您的问题:

Event                Count      Time (sec)     Flops                             --- Global ---  --- Stage ---   Total
                   Max Ratio  Max     Ratio   Max  Ratio  Mess   Avg len Reduct  %T %F %M %L %R  %T %F %M %L %R Mflop/s
------------------------------------------------------------------------------------------------------------------------

--- Event Stage 0: Main Stage

VecScale              42 1.0 1.0709e+00 1.0 2.10e+08 1.0 0.0e+00 0.0e+00 0.0e+00  4 50  0  0  0   4 50  0  0  0   392
MatMult               42 1.0 5.7360e+00 1.1 2.10e+08 1.0 0.0e+00 0.0e+00 0.0e+00 20 50  0  0  0  20 50  0  0  0    73

因此，42 个 VecScale() 操作花费了 1 秒，而 42 个 MatMult() 操作花费了 5.7 秒。在最好的情况下，抑制 VecScale() 操作将使代码速度提高 20%。 for 循环造成的开销甚至更低。我想这就是这个函数不存在的原因。

对于我的计算机性能不佳(VecScale() 为 392Mflops...)，我深表歉意。我很想知道你身上发生了什么!