petsc - 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

这个函数（或任何接近的函数）似乎不在Mat 上运行的函数列表中。因此，对您的问题的简短回答是……不。

如果您经常使用 $y=\frac12 Ax$，一个解决方案是使用MatScale(A,0.5);.

这样的功能有用吗？检查这一点的一种方法是使用-log_summarypetsc 选项来获取一些分析信息。如果您的矩阵是密集的，您会发现在 中花费的时间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

所以 42VecScale()次操作需要 1 秒，而 42 次MatMult()操作需要 5.7 秒。在最好的情况下，抑制VecScale()操作将使代码加速 20%。由于 for 循环的开销甚至更低。我想这就是为什么这个功能不存在的原因。

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