initial values of matrix and vector
Barry Smith
bsmith at mcs.anl.gov
Tue Feb 24 15:35:53 CST 2009
Jed is right. Our intention is to have the vector initialized to
zero at creation time IF we allocate the array space.
Tracking down a not properly initialized vector is the type of bug
that takes forever to find; plus I wanted consistency
between vectors and matrix. Users SHOULD NOT call a MatZeroEntries()
initially on matrices (since it screws up
preallocation) so we should not have them do it on vectors either for
consistency.
Barry
On Feb 24, 2009, at 10:57 AM, Jed Brown wrote:
> On Tue 2009-02-24 06:53, Matthew Knepley wrote:
>> On Mon, Feb 23, 2009 at 8:26 PM, Barry Smith <bsmith at mcs.anl.gov>
>> wrote:
>>
>>
>> Vectors have all entries equal to zero. Dense matrices have
>> all entries
>> equal to zero. Sparse matrices have no entries (logically this
>> is the same
>> as entries equal to zero).
>>
>>
>> We do not automatically zero the entires upon allocation. You have
>> to call
>> VecZeroEntries() or VecSet()
>> to initialize the vector.
>
> That's what I would have expected, but Barry is correct, at least for
> the usual types. See, for instance
>
> PetscErrorCode PETSCVEC_DLLEXPORT VecCreate_Seq(Vec V)
> {
> Vec_Seq *s;
> PetscScalar *array;
> PetscErrorCode ierr;
> PetscInt n = PetscMax(V->map->n,V->map->N);
> PetscMPIInt size;
>
> PetscFunctionBegin;
> ierr = MPI_Comm_size(((PetscObject)V)->comm,&size);CHKERRQ(ierr);
> if (size > 1) {
> SETERRQ(PETSC_ERR_ARG_WRONG,"Cannot create VECSEQ on more than
> one process");
> }
> ierr = PetscMalloc(n*sizeof(PetscScalar),&array);CHKERRQ(ierr);
> ierr = PetscLogObjectMemory(V, n*sizeof(PetscScalar));CHKERRQ(ierr);
> ierr = PetscMemzero(array,n*sizeof(PetscScalar));CHKERRQ(ierr);
> ierr = VecCreate_Seq_Private(V,array);CHKERRQ(ierr);
> s = (Vec_Seq*)V->data;
> s->array_allocated = array;
> PetscFunctionReturn(0);
> }
>
>
>
> In well-structured code, the cost of zeroing the vector is tiny and
> the
> effort to track down uninitialized bugs is significant enough that
> this
> seems like a sane default.
>
> Jed
More information about the petsc-users
mailing list