Line data    Source code

       1             : /*
       2             :    - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
       3             :    SLEPc - Scalable Library for Eigenvalue Problem Computations
       4             :    Copyright (c) 2002-, Universitat Politecnica de Valencia, Spain
       5             : 
       6             :    This file is part of SLEPc.
       7             :    SLEPc is distributed under a 2-clause BSD license (see LICENSE).
       8             :    - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
       9             : */
      10             : 
      11             : #include <petsc/private/vecimpl.h>
      12             : 
      13             : #if defined(__WITH_MPI__)
      14             : #define __SUF__(A) A##_MPI
      15             : #else
      16             : #define __SUF__(A) A##_Seq
      17             : #endif
      18             : #define __QUOTEME_(x) #x
      19             : #define __QUOTEME(x) __QUOTEME_(x)
      20             : #define __SUF_C__(A) __QUOTEME(__SUF__(A))
      21             : 
      22       21764 : static PetscErrorCode __SUF__(VecDot_Comp)(Vec a,Vec b,PetscScalar *z)
      23             : {
      24       21764 :   PetscScalar    sum = 0.0,work;
      25       21764 :   PetscInt       i;
      26       21764 :   Vec_Comp       *as = (Vec_Comp*)a->data,*bs = (Vec_Comp*)b->data;
      27             : 
      28       21764 :   PetscFunctionBegin;
      29       21764 :   SlepcValidVecComp(a,1);
      30       21764 :   SlepcValidVecComp(b,2);
      31       21764 :   if (as->x[0]->ops->dot_local) {
      32       43907 :     for (i=0,sum=0.0;i<as->n->n;i++) {
      33       22143 :       PetscUseTypeMethod(as->x[i],dot_local,bs->x[i],&work);
      34       22143 :       sum += work;
      35             :     }
      36             : #if defined(__WITH_MPI__)
      37       21764 :     work = sum;
      38       65292 :     PetscCallMPI(MPIU_Allreduce(&work,&sum,1,MPIU_SCALAR,MPIU_SUM,PetscObjectComm((PetscObject)a)));
      39             : #endif
      40             :   } else {
      41           0 :     for (i=0,sum=0.0;i<as->n->n;i++) {
      42           0 :       PetscCall(VecDot(as->x[i],bs->x[i],&work));
      43           0 :       sum += work;
      44             :     }
      45             :   }
      46       21764 :   *z = sum;
      47       21764 :   PetscFunctionReturn(PETSC_SUCCESS);
      48             : }
      49             : 
      50       16701 : static PetscErrorCode __SUF__(VecMDot_Comp)(Vec a,PetscInt n,const Vec b[],PetscScalar *z)
      51             : {
      52       16701 :   PetscScalar    *work,*work0,*r;
      53       16701 :   Vec_Comp       *as = (Vec_Comp*)a->data;
      54       16701 :   Vec            *bx;
      55       16701 :   PetscInt       i,j;
      56             : 
      57       16701 :   PetscFunctionBegin;
      58       16701 :   SlepcValidVecComp(a,1);
      59       63442 :   for (i=0;i<n;i++) SlepcValidVecComp(b[i],3);
      60             : 
      61       16701 :   if (as->n->n == 0) {
      62           0 :     *z = 0;
      63           0 :     PetscFunctionReturn(PETSC_SUCCESS);
      64             :   }
      65             : 
      66       16701 :   PetscCall(PetscMalloc2(n,&work0,n,&bx));
      67             : 
      68             : #if defined(__WITH_MPI__)
      69       16701 :   if (as->x[0]->ops->mdot_local) {
      70       16701 :     r = work0; work = z;
      71             :   } else
      72             : #endif
      73             :   {
      74           0 :     r = z; work = work0;
      75             :   }
      76             : 
      77             :   /* z[i] <- a.x' * b[i].x */
      78       63442 :   for (i=0;i<n;i++) r[i] = 0.0;
      79       33660 :   for (j=0;j<as->n->n;j++) {
      80       64216 :     for (i=0;i<n;i++) bx[i] = ((Vec_Comp*)b[i]->data)->x[j];
      81       16959 :     if (as->x[0]->ops->mdot_local) PetscUseTypeMethod(as->x[j],mdot_local,n,bx,work);
      82           0 :     else PetscCall(VecMDot(as->x[j],n,bx,work));
      83       64216 :     for (i=0;i<n;i++) r[i] += work[i];
      84             :   }
      85             : 
      86             :   /* If def(__WITH_MPI__) and exists mdot_local */
      87             : #if defined(__WITH_MPI__)
      88       16701 :   if (as->x[0]->ops->mdot_local) {
      89             :     /* z[i] <- Allreduce(work[i]) */
      90       50103 :     PetscCallMPI(MPIU_Allreduce(r,z,n,MPIU_SCALAR,MPIU_SUM,PetscObjectComm((PetscObject)a)));
      91             :   }
      92             : #endif
      93             : 
      94       16701 :   PetscCall(PetscFree2(work0,bx));
      95       16701 :   PetscFunctionReturn(PETSC_SUCCESS);
      96             : }
      97             : 
      98           3 : static PetscErrorCode __SUF__(VecTDot_Comp)(Vec a,Vec b,PetscScalar *z)
      99             : {
     100           3 :   PetscScalar    sum = 0.0,work;
     101           3 :   PetscInt       i;
     102           3 :   Vec_Comp       *as = (Vec_Comp*)a->data,*bs = (Vec_Comp*)b->data;
     103             : 
     104           3 :   PetscFunctionBegin;
     105           3 :   SlepcValidVecComp(a,1);
     106           3 :   SlepcValidVecComp(b,2);
     107           3 :   if (as->x[0]->ops->tdot_local) {
     108           9 :     for (i=0,sum=0.0;i<as->n->n;i++) {
     109           6 :       PetscUseTypeMethod(as->x[i],tdot_local,bs->x[i],&work);
     110           6 :       sum += work;
     111             :     }
     112             : #if defined(__WITH_MPI__)
     113           3 :     work = sum;
     114           9 :     PetscCallMPI(MPIU_Allreduce(&work,&sum,1,MPIU_SCALAR,MPIU_SUM,PetscObjectComm((PetscObject)a)));
     115             : #endif
     116             :   } else {
     117           0 :     for (i=0,sum=0.0;i<as->n->n;i++) {
     118           0 :       PetscCall(VecTDot(as->x[i],bs->x[i],&work));
     119           0 :       sum += work;
     120             :     }
     121             :   }
     122           3 :   *z = sum;
     123           3 :   PetscFunctionReturn(PETSC_SUCCESS);
     124             : }
     125             : 
     126           3 : static PetscErrorCode __SUF__(VecMTDot_Comp)(Vec a,PetscInt n,const Vec b[],PetscScalar *z)
     127             : {
     128           3 :   PetscScalar    *work,*work0,*r;
     129           3 :   Vec_Comp       *as = (Vec_Comp*)a->data;
     130           3 :   Vec            *bx;
     131           3 :   PetscInt       i,j;
     132             : 
     133           3 :   PetscFunctionBegin;
     134           3 :   SlepcValidVecComp(a,1);
     135           9 :   for (i=0;i<n;i++) SlepcValidVecComp(b[i],3);
     136             : 
     137           3 :   if (as->n->n == 0) {
     138           0 :     *z = 0;
     139           0 :     PetscFunctionReturn(PETSC_SUCCESS);
     140             :   }
     141             : 
     142           3 :   PetscCall(PetscMalloc2(n,&work0,n,&bx));
     143             : 
     144             : #if defined(__WITH_MPI__)
     145           3 :   if (as->x[0]->ops->mtdot_local) {
     146           3 :     r = work0; work = z;
     147             :   } else
     148             : #endif
     149             :   {
     150           0 :     r = z; work = work0;
     151             :   }
     152             : 
     153             :   /* z[i] <- a.x' * b[i].x */
     154           9 :   for (i=0;i<n;i++) r[i] = 0.0;
     155           9 :   for (j=0;j<as->n->n;j++) {
     156          18 :     for (i=0;i<n;i++) bx[i] = ((Vec_Comp*)b[i]->data)->x[j];
     157           6 :     if (as->x[0]->ops->mtdot_local) PetscUseTypeMethod(as->x[j],mtdot_local,n,bx,work);
     158           0 :     else PetscCall(VecMTDot(as->x[j],n,bx,work));
     159          18 :     for (i=0;i<n;i++) r[i] += work[i];
     160             :   }
     161             : 
     162             :   /* If def(__WITH_MPI__) and exists mtdot_local */
     163             : #if defined(__WITH_MPI__)
     164           3 :   if (as->x[0]->ops->mtdot_local) {
     165             :     /* z[i] <- Allreduce(work[i]) */
     166           9 :     PetscCallMPI(MPIU_Allreduce(r,z,n,MPIU_SCALAR,MPIU_SUM,PetscObjectComm((PetscObject)a)));
     167             :   }
     168             : #endif
     169             : 
     170           3 :   PetscCall(PetscFree2(work0,bx));
     171           3 :   PetscFunctionReturn(PETSC_SUCCESS);
     172             : }
     173             : 
     174       17872 : static PetscErrorCode __SUF__(VecNorm_Comp)(Vec a,NormType t,PetscReal *norm)
     175             : {
     176       17872 :   PetscReal      work[3],s=0.0;
     177       17872 :   Vec_Comp       *as = (Vec_Comp*)a->data;
     178       17872 :   PetscInt       i;
     179             : 
     180       17872 :   PetscFunctionBegin;
     181       17872 :   SlepcValidVecComp(a,1);
     182             :   /* Initialize norm */
     183       17872 :   switch (t) {
     184          12 :     case NORM_1: case NORM_INFINITY: *norm = 0.0; break;
     185       17857 :     case NORM_2: case NORM_FROBENIUS: *norm = 1.0; s = 0.0; break;
     186           3 :     case NORM_1_AND_2: norm[0] = 0.0; norm[1] = 1.0; s = 0.0; break;
     187             :   }
     188       36067 :   for (i=0;i<as->n->n;i++) {
     189       18195 :     if (as->x[0]->ops->norm_local) PetscUseTypeMethod(as->x[i],norm_local,t,work);
     190           0 :     else PetscCall(VecNorm(as->x[i],t,work));
     191             :     /* norm+= work */
     192       18195 :     switch (t) {
     193           6 :       case NORM_1: *norm+= *work; break;
     194       18165 :       case NORM_2: case NORM_FROBENIUS: AddNorm2(norm,&s,*work); break;
     195           6 :       case NORM_1_AND_2: norm[0]+= work[0]; AddNorm2(&norm[1],&s,work[1]); break;
     196          24 :       case NORM_INFINITY: *norm = PetscMax(*norm,*work); break;
     197             :     }
     198       18195 :   }
     199             : 
     200             :   /* If def(__WITH_MPI__) and exists norm_local */
     201             : #if defined(__WITH_MPI__)
     202       17872 :   if (as->x[0]->ops->norm_local) {
     203       17872 :     PetscReal work0[3];
     204             :     /* norm <- Allreduce(work) */
     205       17872 :     switch (t) {
     206           3 :     case NORM_1:
     207           3 :       work[0] = *norm;
     208           9 :       PetscCallMPI(MPIU_Allreduce(work,norm,1,MPIU_REAL,MPIU_SUM,PetscObjectComm((PetscObject)a)));
     209             :       break;
     210       17857 :     case NORM_2: case NORM_FROBENIUS:
     211       17857 :       work[0] = *norm; work[1] = s;
     212       53571 :       PetscCallMPI(MPIU_Allreduce(work,work0,1,MPIU_NORM2,MPIU_NORM2_SUM,PetscObjectComm((PetscObject)a)));
     213       17857 :       *norm = GetNorm2(work0[0],work0[1]);
     214       17857 :       break;
     215           3 :     case NORM_1_AND_2:
     216           3 :       work[0] = norm[0]; work[1] = norm[1]; work[2] = s;
     217           9 :       PetscCallMPI(MPIU_Allreduce(work,work0,1,MPIU_NORM1_AND_2,MPIU_NORM2_SUM,PetscObjectComm((PetscObject)a)));
     218           3 :       norm[0] = work0[0];
     219           3 :       norm[1] = GetNorm2(work0[1],work0[2]);
     220           3 :       break;
     221           9 :     case NORM_INFINITY:
     222           9 :       work[0] = *norm;
     223          27 :       PetscCallMPI(MPIU_Allreduce(work,norm,1,MPIU_REAL,MPIU_MAX,PetscObjectComm((PetscObject)a)));
     224             :       break;
     225             :     }
     226       17872 :   }
     227             : #else
     228             :   /* Norm correction */
     229           0 :   switch (t) {
     230           0 :     case NORM_2: case NORM_FROBENIUS: *norm = GetNorm2(*norm,s); break;
     231           0 :     case NORM_1_AND_2: norm[1] = GetNorm2(norm[1],s); break;
     232           0 :     default: ;
     233             :   }
     234             : #endif
     235       17872 :   PetscFunctionReturn(PETSC_SUCCESS);
     236             : }
     237             : 
     238         183 : PetscErrorCode __SUF__(VecDotNorm2_Comp)(Vec v,Vec w,PetscScalar *dp,PetscScalar *nm)
     239             : {
     240         183 :   PetscScalar       dp0=0.0,nm0=0.0,dp1=0.0,nm1=0.0;
     241         183 :   const PetscScalar *vx,*wx;
     242         183 :   Vec_Comp          *vs = (Vec_Comp*)v->data,*ws = (Vec_Comp*)w->data;
     243         183 :   PetscInt          i,n;
     244         183 :   PetscBool         t0,t1;
     245             : #if defined(__WITH_MPI__)
     246          90 :   PetscScalar       work[4];
     247             : #endif
     248             : 
     249         183 :   PetscFunctionBegin;
     250             :   /* Compute recursively the local part */
     251         183 :   PetscCall(PetscObjectTypeCompare((PetscObject)v,VECCOMP,&t0));
     252         183 :   PetscCall(PetscObjectTypeCompare((PetscObject)w,VECCOMP,&t1));
     253         183 :   if (t0 && t1) {
     254          90 :     SlepcValidVecComp(v,1);
     255          90 :     SlepcValidVecComp(w,2);
     256         183 :     for (i=0;i<vs->n->n;i++) {
     257          93 :       PetscCall(VecDotNorm2_Comp_Seq(vs->x[i],ws->x[i],&dp1,&nm1));
     258          93 :       dp0 += dp1;
     259          93 :       nm0 += nm1;
     260             :     }
     261          93 :   } else if (!t0 && !t1) {
     262          93 :     PetscCall(VecGetLocalSize(v,&n));
     263          93 :     PetscCall(VecGetArrayRead(v,&vx));
     264          93 :     PetscCall(VecGetArrayRead(w,&wx));
     265        9679 :     for (i=0;i<n;i++) {
     266        9586 :       dp0 += vx[i]*PetscConj(wx[i]);
     267        9586 :       nm0 += wx[i]*PetscConj(wx[i]);
     268             :     }
     269          93 :     PetscCall(VecRestoreArrayRead(v,&vx));
     270          93 :     PetscCall(VecRestoreArrayRead(w,&wx));
     271           0 :   } else SETERRQ(PetscObjectComm((PetscObject)v),PETSC_ERR_ARG_INCOMP,"Incompatible vector types");
     272             : 
     273             : #if defined(__WITH_MPI__)
     274             :     /* [dp, nm] <- Allreduce([dp0, nm0]) */
     275          90 :     work[0] = dp0; work[1] = nm0;
     276         270 :     PetscCallMPI(MPIU_Allreduce(work,&work[2],2,MPIU_SCALAR,MPIU_SUM,PetscObjectComm((PetscObject)v)));
     277          90 :     *dp = work[2]; *nm = work[3];
     278             : #else
     279          93 :     *dp = dp0; *nm = nm0;
     280             : #endif
     281         183 :   PetscFunctionReturn(PETSC_SUCCESS);
     282             : }
     283             : 
     284             : #undef __SUF__
     285             : #undef __QUOTEME
     286             : #undef __SUF_C__