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()上花费的时间大得多。只有在处理了一个sparce矩阵时，这个问题才有意义，每行有几个非空项。

下面是一个代码进行测试，使用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循环的开销甚至比这更低。我想这就是为什么这个函数不存在的原因。

我很抱歉我的电脑性能不佳（392Mflops for VecScale（）...）。我很想知道你的情况。